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EVALUATION OF SUBGRID-SCALE TURBULENCE MODELS 
USING A FULLY SIMULATED TURBULENT FLOW 

Abstract 


Numerous models have been proposed for approximating the subgrid- 
scale Reynolds stresses m numerical simulations of turbulent fluid flow. 
Until now, the only way to verify such approximations has been to observe 
the resulting behavior of the large-scale flow. If the entire turbulent 
flow field were known, it would be possible to make direct comparisons 
between the exact Reynolds stresses and a given model. We have calcula- 
ted an "exact’ 1 turbulent flow field on a three-dimensional grid with 64 
points on a side. The flow simulates grid-generated turbulence from x?ind 
tunnel experiments. In this simulation, the grid spacing is small enough 
to include essentially all of the viscous energy dissipation and the box 
is large enough to contain the largest eddy m the flow. The method is 
limited to low-turbulence Reynolds numbers, in our case R^ = 36.6. 

In order to complete the calculation using a reasonable amount of 
computer time with reasonable accuracy, we developed a third-order time- 
integration scheme which runs at about the same speed as a simple first- 
order scheme. It obtains this accuracy by saving the velocity field and 
its first-time derivative at each time step. Fourth-order accurate space- 
differencmg is used. 

The results of this simulation were treated as an experimental reali- 
zation of physical turbulence. We then superimposed an 8x8x8 coarse 
mesh over the original fine mesh and defined a filtered velocity field 


u^Cx) as the local spatial average of u . From these we defined the 
subgrid-scale velocity field u^ by u = u + u^. The filtering process 
gives rise to three terms in the Navier-Stokes equations. These are the 


Reynolds stress, u'u', the Leonard term, u u - _u u , and the cross 
— — -f — ij i] l J 

term, u'u + u u'. We demonstrated that the cross term is non-zero and 
3 i J 

is, m fact, dissipative; we also developed a model for it. The Leonard 


term and the cross term can be combined into a single term which can be 


modeled by (Vu )*(Vu ). This reduces, in one dimension, to a quadratic 
a J 

artificial viscosity frequently used in compressible flow calculations 
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The relationship between the filtering process and artificial viscosity 
is shown. 

Finally we calculated each of the above terms within each cell on 
the coarse mesh, and we attempted to model them using the f-i-ltered veloc- 
ity field. For each model we calculated the correlation between the model 
and its "exact" value. We found the correlation between the Leonard and 
cross terms and their models to be excellent, around ninety percent. The 
correlation between model and experiment for the Reynolds stress is not 
as good, but we did achieve about seventy percent correlation between the 
dissipation produced by the Reynolds stress and its model. We found no 
model that is significantly better than the standard Smagor insky model. 

We found that models using the subgr id-scale turbulent kinetic energy are 
no better than Smagor insky ’ s , even when we had the exact subgrid-scale 
kinetic energy to work with. All of these conclusions must be qualified 
by stating that we were working at very low turbulent Reynolds numbers, 
and the results cannot necessarily be extrapolated to higher Reynolds num- 
bers. 
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Chapter I 
INTRODUCTION 


We begin with a brief discussion of the general approach to the 
numerical simulation of turbulent flows. It is generally not 'possible to 
calculate a turbulent flow in complete detail, because the range of length 
scales involved is so large that the amount of data that would have to be 
handled is orders of magnitude greater than the capacity of any existing 
or projected computer. For this reason, the traditional approaches to 
such problems have been based on Reynolds’ original idea of averaging the 
Navier-Stokes equations over an ensemble of identical flows or an approp- 
riate interval of time or space. One then has equations for an averaged 
velocity field u{x,t), where the overbar denotes averaging according to 
whatever definition is employed. If we then define a fluctuating velocity 
component by u(x,t) = u(x,t) + u’(x,t), the averaged equations can be 
written (for an incompressible flow with constant viscosity) 


3u 

i 

at 


a — 

T U U. = 

ax. i j 


= _-i£ 


+ vV u 
3x. i 

l 


1ST 


ax 


= 0 


( 1 . 1 ) 


( 1 . 2 ) 


R = u’u’ , (1.3) 

ij id 

since u = u. and u’ = 0 as a consequence of the definition of aver- 
ill 

aging. To formally close this system of equations it is necessary to find 
an expression for R (the Reynolds stress) m terms of u^. 

The search for such expressions ( closure models) has been a major 
direction m turbulence work for many decades. Prior to the advent of 
computers, only very simple models could be used, i.e., those which yielded 
equations to which one could obtain solutions analytically. These models, 
which assume that the Reynolds stresses, like viscous stresses, are propor- 
tional to the local strain rate of the fluid, give acceptable results 
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provided the range of flows they are required to predict is not too 
large These models ( eddy viscosity models ) can be adjusted for partic- 
ular flows to produce excellent results. The problem is that this type 
of model must be adjusted to experimental dat-a and can therefore be used 
only in an interpolative manner 

More advanced models have been proposed since the introduction of 
large digital computers. These include models in which the eddy vis- 
cosity is a function of space and/or time ( turbulence kinetic energy and 
two- equation models) and still more advanced models which utilize par- 
tial differential equations for the Reynolds stresses are currently being 
developed. It is too early to predict the success of these approaches, 
but there is reason to believe that their range of application will be 
limited. By this we mean that the model parameters will probably need 
to be adjusted for each type of flow. 

To see why this might be the case and what might be done about it, 
consider the overall nature of a turbulent flow field. The range in 
length scales between the largest and smallest turbulent structures is 
many orders of magnitude in. most flows of interest The largest turbu- 
lent structures draw energy from the mean flow. This energy is thought 
to cascade through an intermediate range of eddy sizes to eventually 
reach the smallest turbulent eddies. The smallest eddies then dissipate 
the kinetic energy to internal energy by viscous effects. There is rea- 
son to believe that the largest structures m a turbulent flow are much 
more dependent on the origin of the turbulence, i.e., the type of flow 
under consideration, than either the intermediate or small-scale struc- 
tures. This could explain the failure of any single model to predict a 


wide variety of flows when the large-scale turbulence is included in 

u'u'. It appears more likely that a universal model might be found if 
i J 

only the intermediate and small-scale turbulence is modeled. This ap- 


proach ( large eddy simulation ) requires that the large-scale turbulence 
be calculated explicitly and that the small scales (the "subgrid" scale 


turbulence) be modeled. With the recent advances in high-speed computing 


machinery, this approach has become increasingly practical m the last 


several years (Hirt, 1969; Deardorff, 1970; Fox and Lilly, 1972). 
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In large eddy simulation, one averages the Navier-Stokes equations 
over a small spatial region m order to remove the small-scale fluctua- 
tions. The resulting equation for the large-scale field contains a term 
similar to, but more complicated than, the Reynolds stress of Eqs. 

(1 1) and (1.3), and this term (the subgrid scale Reynolds stress) must 
be modeled. 

Several models for the subgrid scale R have been proposed. The 
problem has been how to verify a proposed model. The best that could be 
done until now was to compare the evolution of the large-scale structures 
in a computation to those m an experiment. This will not reveal whether 
or not the actual subgrid scale Reynolds stress is being accurately mod- 
eled, but only whether or not the subgrid scale Reynolds stress and the 
model have the same net effect on the large-scale motions for the particu- 
lar type of flow m question. In addition, virtually all models contain 
at least one adjustable constant which must be set by some ad hoc assump- 
tion or by adjusting the constant to fit some important aspect of an ex- 
periment. On the other hand, if there wer'e a physical experiment which 
measured everything of interest m a turbulent flow field, from the larg- 
est turbulent structure to the smallest eddy, it would then be possible 
to compute the subgrid scale R exactly, and then compare its value at 
each point m space to the prediction of a model. Unfortunately, there 
is no laboratory experiment capable of such measurements. But if we could 
compute the evolution of a turbulent flow field on a sufficiently fine 
grid (fine enough to include all of the turbulent structures) numerically, 
then we would have all the information necessary to make direct compari- 
sons between measured values of the subgrid R and the model predictions. 

The objectives of this work were to accurately calculate a three- 
dimensional turbulent flow field on a fine grid by directly integrating 
the Navier-Stokes equations using no approximations with respect to the 
structure of the turbulence, i.e., without having to average the equations, 
and then to use the results of that calculation to examine subgrid scale 
models on a coarse mesh overlayed on the original fine mesh. Practical 
limitations require that this be done at a low Reynolds number, since at 
high Reynolds numbers the difference m scale between the largest and 
smallest eddies is so great that computer simulation is impractical. Also, 
the flow should be as physically simple as possible (see below) 
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In Chapter II we describe the numerical integration method which was 
used to compute the flow field on the fine mesh. A third-order time 
scheme is introduced which has not previously been used m this applica- 
tion Because this is the. first attempt at model verification, we have 
chosen the simplest possible turbulent flow field for our calculation, 
homogeneous isotropic turbulence This avoids any problems of anisotropy, 
but restricts the results to problems in which the subgrid scale turbu- 
lence can be assumed isotropic In Chapter III we describe m detail how 
the mam calculation was performed, including some programming techniques 
we employed to greatly speed up the calculation In Chapter IV we discuss 
the general problem of modeling subgrid scale turbulence with emphasis on 
methods which could be verified by our calculations. In Chapter V we 
demonstrate that the results of our main calculation do m fact have the 
properties of real turbulent flow. We then make the comparisons of the 
numerically calculated values of R and the predictions of various 
models and show that, although the models currently used are not as accu- 
rate as one would like, they are difficult to improve upon in a simple 
manner. 


4 



Chapter II 


THE NUMERICAL METHOD 


2.1 General 

The finite-difference scheme chosen for the main calculation is 
fourth-order accurate m space and third-order accurate in time. The 
rationale for this choice is worth a brief discussion. 

In many large computer simulations a major effort is made to keep 
the problem small enough to be contained entirely m the mam memory of 
the computer. This is done to avoid the use of the disk memory and its 
relatively slow rate of data transfer. The problem we shall try to solve 
has three velocity components at each of 262,144 grid points, for a total 
of 786,432 words necessary to specify the velocity field at one time step. 
This number alone exceeds the roughly 400,000 words of memory available 
xn our CDC 7600 large-core memory. Since we are forced to utilize disk 
memory, waiting times for the completion of data transfer to and from 
disk become a major consideration. Large amounts of data must be trans- 
ferred from disk to mam memory, processed m mam memory, then transferred 
back to disk. If the processing time is too short, the data transfer time 
will determine the running time of the problem. In our case, using 
fourth-order differencing m space, we found that only about five percent 
of the total running time was spent m waiting for data transfer to be 
completed. This means that the data were processed slightly faster than 
it could be transferred, even though a double-buffering scheme was used 
Had we used second-order space differencing, we would not have gained any- 
thing m running time, since the data-transf er time would still be the 
same Reducing the processing time would have simply increased the per- 
centage of wait time. This means that we have used fourth-order space 
differencing with no increase m computer charges with respect to second- 
order differencing, i.e., increased accuracy is obtained at essentially 
no cost. We emphasize that this choice is not made for accuracy reasons, 
but is a result of being forced to use the disk memory. 
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Having settled on fourth-order accuracy in space, we must next decide 
how to handle the time differencing. In any numerical method, common 
sense dictates that the truncation error due to the time differencing 
should be roughly the same as the truncation error due to the space dif- 
ferencing This can be done even with first-order time differencing if 
the time step used is sufficiently small. The criterion for choosing a 
time-differencing scheme now becomes cost, i.e. , computer running time. 

If a second-order time scheme were to take twice the computing time per 
time step as a first-order scheme but allowed us to increase the time step 
by more than a factor of two, the total running time would be reduced and 
the second-order scheme would be justified. 

You generally get what you pay for with numerical methods. If greater 
accuracy is desired, you must pay for it. However, there can sometimes be 
more than one method of payment The most common method of payment is in 
increased running time. For example, the simplest two-step, second-order 
schemes (Roache, 1972) obtain second-accuracy by essentially performing a 
simple first-order scheme twice, thus doubling the running time. Another 
method of payment can be in increased storage requirements. The leap-frog 
scheme performs essentially the same calculations as a first-order scheme, 
but it obtains second-order accuracy by saving an extra time step m the 
calculation, hence doubling the storage requirements. On most present- 
day computing systems, the extra charge for doubling the storage require- 
ment is small m comparison to the savings resulting from halving the CPU 
time. On this basis the second-order accuracy of the leap-frog method 
appears to be obtained almost for free. Another way to look at it is that 
the user is usually charged for most of the available storage, and if he 
does not use it he is short-changing himself. 

With the above in mind, we note that since our problem requires use 
of the disk the storage available is, for all practical purposes, unlimi- 
ted. As noted above, the leap-frog scheme obtains its second-order accu- 
racy by saving the velocity field at an extra time step. We have developed 
a time-diff erencmg schdme that obtains third-order accuracy by saving the 
first time derivative of the velocity field, as well as the velocity field 
itself. With this method the running time per time step is essentially 
unchanged from first-order methods, of from the leap-frog method, but the 
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time step can be increased while maintaining the same truncation error. 
This reduces the total running time of the problem without reducing the 
accuracy. 


2.2 The Tune-Differencing Scheme 

-k 

We now develop the third-order method that we will use. It is essen- 
tially a predictor-corrector method with a second-order leap-frog predic- 
tor and a third-order Adams-Moulton corrector. To explain the method we 
deal with the ordinary differential equation 


u = au 
t 


(2 1 ) 


which has the exact solution u = exp (at). Suppose that we are given u 
at times nSt and (n-l)6t to third-order accuracy, i.e., 


u 11 = exp[an6t] + 0(<5t ) 


u 11 1 = exp [a (n-1) St] + O(dt^) 


(2.2a) 

(2.2b) 


and (u ) at time n<St and (n-l)6t to second-order accuracy, i.e. , 


(u ) n = a exp[an<$t] + 0(6t^) 
t 


(u )" 1 = a exp [a (n-l)<5t] + 0(6t^) 


(2.3a 

(2.3b) 


Thus, (u ) n and u 11 represent two numerical approximations to the exact 
solution; they are, respectively, the predicted and corrected values. Let 
us first approximate the solution to equation (2.1) at time (n+l)6t 
using the standard leap-frog method. 


, * n+1 n-1 , . * N n 

(u ) = u + 26t(u ) 


(2.4) 


Using a Taylor series expansion to get exp (a(n+l)<St) in terms of 
exp(andt) and introducing the notation y = aSt, we have 
* 

For a general discussion of methods of this type, see Appendix A 

1 



exp[ct(n±l)St] 


‘ ( x± 


2 3 4 

v + — + 2— + X_ 
Y 2 “ 6 24 




exp[an6t] (2.5) 


from which it is easy to show that 


* n-J-1 v 4 

(u ) = exp[a(n+l)6t] - exp [anSt] + 0(<5t ) (2.6) 


( *,n+l * n+1 , * n+1 

Hence, (u ) and (u ) = a(u ) are accurate to second-order. 

* n+1 * n * n— 1 

Now that we have (u ) , (u ) , and (u ) to second-order, we 

^ Tl ^ ^ * C n 

can evaluate (u ) to first-order and (u ) to zeroth order. 

tt ttt 

f *\ n+1 , * n-1 

(u h - (u >t 2 2 

= a exp(andt) + 0(<5t ) (2.7) 

2<$t 

(u ) - 2(u ) + (u ) 

E — = a J exp ( an( s t ) + 0 (6t) (2.8) 


(u ) 


tt 


(u ) 


n 

ttt 


fit 


Next we evaluate the corrected value u n3 "^~ to third-order accuracy by 
using the expansion (2.5) m the form 


n+1 n * n St 2 , * n St 3 , * n 

u - u + «t(u ) t +— (u ) tfc +-g- (u ) 


ttt 


(2.9) 


which is identical to 


■i n+1 = [exp (an<St) + 0(6t 4 )] + 5t[a exp(anfit) + 0(6t 3 )] (2.10) 

6t 2 2 2 « 

+ ~ 2 ~ exp(anfit) + 0(6t )] H [a exp(anSt) + 0(<$t)] 


Comparing (2.5) and (2.10), we see that 


u n+1 = exp [a(n+l)<$t] + 0(6t^) 


Summarizing the method, the predictor step is leap frog 


(uV +1 - u n_1 + 2«t(u*)" 


(2 11 ) 


and the corrector step m simplified form is 
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(2 12 ) 


n+1 

u 


n 


u + ~ <5t(u + ~ 6t(u 



n-1 

t 


which will he recognized as an implicit Ad a ms -Moult on method. It is, m 
fact, possible to produce a fourth-order method m this way, but the extra 
advantage is minimal. 


2.3 Numerical Stability of the Third-Order Method 

Following standard Von Neumann analysis (Richtmyer, 1967), we seek 
solutions to equations (2.11) and (2.12) of the form 


u n = A n u (2.13a) 

* n n * 

(u ) n = Au (2.13b) 


where A is m general a complex constant. The numerical method will be 
stable if we can ensure that j A | <_ 1 for Re(cx) 0, First, we replace 
(u by a(u ) n m Eq. (2.11) to get 

n 1 / nt 1 n / n > 

u = (u ) - 2y(u ) (2.14) 

where again y = a<5t. Changing the index, we have 

n . n+2 0 . *, n+1 , 9 . r 

u — (u ) — 2y (u ) (2 15a) 


and 


n+1 

u 


, *•, n+3 „ . * n+2 

(u ) - 2y(u ) 


(2 15b) 


-fj 

Substituting (2.15a) and (2.15b) into (2 12) and again replacing (u ) 

* n Z 

by a(u ) yields 


* n+'l * n+2 

(u ) n+J - 2y(u ) n+Z 


*, n+2 


* n+1 


(u ■')““■ - 2y(u )“ Ti + f y(u ) 

, 5 . * n+1 1 , * n-1 

+ 12 Y U ' - ^ Y(u ) 


n 


(2 16) 


9 



yields a 


Substituting (2.13b) into (2.16) and multiplying by A n+ ^ 
quartic equation for A: 

L ■j IQ o 9 v 

A* ^ (l+2y)A J + g yA - fY A + YJ = 0 (2.17) 

The stability of the numerical solution is now determined by finding 
the maximum value of |y | = |a|dt with Re(a) < 0 for which all four 
roots of Eq. (2.17) have magnitudes less than 1. This ensures that the 
solutions given by Eq. (2 13) will not increase exponentially with in- 
creasing n m cases in which they should not. The resulting region of 
stability is shown, m Fig. 2.1. The curve shown m Fig. 2.1 is the curve 
of neutral stability on which the magnitude of the largest A is 1. 

This curve was found by computing the roots of Eq. (2.17) for fixed y^ 
noting the value of |y^_ | for which one of the roots becomes unity. 

Other methods are available, but the simplicity of this calculation does 
not warrant their use. 

2.4 Accuracy of the Third-Order Method 

We again consider the ordinary differential equation 

u fc = ctu (2.18) 


which has the exact solution 

u = exp (at) = exp(yn) (2.19) 

One of the roots of Eq. (2.17) will be an accurate approximation to 
exp(y). The others are parasitic or computational roots. It turns out 
that the three parasitic roots are highly damped (y^-negative and large 
m comparison to the non-parasitic root) , so that the solution obtained 
will be 


u (n<$t) = A^ 


( 2 . 20 ) 
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Fig. 2.1 



where A is the desired root of Eq (2.17). Let y = Re(y) and 

JL 

y = Im(y). The accuracy of the solution is determined by how well A^ 
approximates exp(y). We separate the numerical error into the phase 
and amplitide errors. The phase error is given by 

Ay = (y^ “ Im[Jln(A^) ] } (2.21) 

and the amplitude error is given by 


Ay t = (y r - Re[£n(A 2 )]} 


( 2 . 22 ) 


We have computed A^ for the leap-frog and third-order schemes. In 
Table 2.1 we list some values of y^ and the error m the computed value 
of y , Ay^, for the leap-frog and third-order schemes with y r = 0. In 
Table 2.2 we list some values of y r and the error in the computed value 
of y , Ay , for the two schemes with y =0. The range of y^ and y 

1 IT 1 1 L 

listed in the tables covers the range of interest in our mam calculation. 

A simple linear equation which is sometimes used as a model for test- 
ing numerical approximations to the Navier-Stokes equations is the convec- 
ted diffusion equation: 


9u . 8u 
8t + c 3x 


v 



Assuming a solution to this equation of the form u(x,t) 
have 


lkx 

u(t) e , we 


u. 


= (-ick - vk )u 


(2.23) 


The values for mean velocity, viscosity, spatial increment, and time step 

which correspond to the mam calculation described m Chapter III are 

2 

c = 5.5 cm/sec, v = 0.14 cm /sec, 6x = 20/64 cm, 6t = 0.0073 sec, and 

0 < k < ir/6x Using these values, we have calculated y = -ckdt and 

2 1 
y^ = — vk 6t for use m Tables 2.1 and 2.2. 

The Ay's listed m the tables are the errors per time step. Take, 

for example, y = 0.201 (this is equivalent to a wave with wavelength 
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Table 2.1 


Phase Error Comparison 




2nd 

Order 

3rd 

Order 

K 

^i 

Ay 

i 

a Yi 

.503 

020 

1 4 x 10 -6 

7.6 x IO" 9 

1.005 

.040 

1.1 x 10“ 5 

2 5 x IO" 8 

1.508 

040 

3.7 x 10“ 5 

1.9 x IO' 7 

2.011 

.080 

8.8 x 10 -5 

8.0 x 10“ 7 

2 513 

.100 

1.7 x lO" 4 

2.4 x 10 -6 

3.016 

.120 

3.0 x 10“ 4 

6.0 x IO -6 

3.519 

.141 

4 8 x 10~ 4 

1.2 x 10~ 5 

4 021 

.161 

7.2 x 10“ 4 

2 5 x 10“ 5 

4.524 

.181 

1.0 x io -3 

4.4 x 10“ 5 

5.027 

.201 

1.4 x IO -3 

7.6 x 10“ 5 

5.529 

.221 

1.8 x 10 -3 

1.2 x IO" 4 

6.032 

.241 

2 4 x 10 -3 

1.8 x 10 -4 

6.535 

261 

3.1 x 10 -3 

2 7 x 10“ 4 

7.037 

.281 

3.8 x 10 -3 

3.8 x 10“ 4 

7 540 

.302 

4.8 x 10 -3 

5.2 x 10“ 4 

8.042 

.322 

6.0 x 10“ 3 

7.6 x IO" 4 

8.545 

.342 
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1.0 x IO" 3 

9 048 

.362 

8 4 x 10~ 3 

1.3 x IO” 3 

9.550 

380 

1 0 x io -2 

1.7 x 10“ 3 

10.053 

402 

1.2 x 10“ 2 

2.1 x 10“ 3 
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Table 2.2 


Amplitude Error Comparison 
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1.8 x 

io“ 4 

3 5 

X 

io“ 5 

2.4 x 

10“ 4 

5.3 

X 

10“ 5 
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4^x in our mam calculation). For y = 0.201, the phase error per 
time step is 18 times larger for the leap-frog scheme than for the third- 
order scheme. This means that if both schemes used the same time step the 
accumulated error at a given point m time would be 18 times greater using 
leap frog. This is not the whole story The relevant question to ask is 
how much would you have to reduce the time step using the leap-frog scheme 
m order to get the same error as with the third-order scheme? Suppose we 
are using the third-order scheme with a <5t such that y = 0.201; this 
value is typical of the problem that we actually ran. If we reduce the 
time step by a factor of 4.4 and use the leap-frog scheme, the phase error 

_5 

per time step will be 1.9 * 10 , but we will have to run 4.4 times as 

many time steps so that the phase error per original time step will be 
-5 -5 

4.4 x 1.9 x 10 =8.36x10 , which is slightly larger than the 

7.6 x 10 ^ for the third-order scheme. For smaller y the factors are 

larger than 4.4, and for larger y the factors are smaller than 4.4 Our 
conclusion is that for the leap-frog scheme to achieve the same accuracy 
as the third-order scheme the time step would have to be reduced by at 
least a factor of four, thus increasing the running time of the problem 
by nearly a factor of four. 

To further demonstrate the accuracy of the third-order scheme, we pre- 
sent, in Fig. 2.2, the results of a numerical test on the one-dimensional 
wave equation: 


u + cu =0 (2.24) 

u^ was calculated by the use of Fourier transforms so that the only error 
in the numerical solution is due to the time-dif f erencmg scheme, u was 
defined at 64 evenly spaced points and was initially zero, except for the 
triangle shown near the center. Periodic boundary conditions were applied, 
and we set u<St/6x =0.2. In the exact solution, the triangle moves from 
left to right at the constant speed c so that after 1600 time steps the 
triangle should have swept across the grid five times and the exact solu- 
tion is identical to the initial conditions The third-order calculation 
used 3% more computing time than the leap-frog calculation and twice the 
amount of storage. The starting conditions m each case were exact This 
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was done by calculating the Fourier transform of the initial conditions, 
multiplying each Fourier component by the appropriate A^(ick6t) and 
inverting the resulting Fourier field to get the value of u at time St. 

2.5 Starting the Third-Order Method 

Given the initial conditions for a problem at time t = 0, we can- 
not get the solution at time t = 6t by the third-order method. The 
field at t = St must be found by another method. The technique we used 
is to divide the initial time step St into three increments of St/3 
and use a two-step predictor corrector method to obtain u(St) to second- 
order accuracy. We tested this, method by again solving the equation 
u t + cu x =0 as above, except that rather than using the exact numerical 
solution for t = St we used the second-order predictor corrector scheme 
to get the second time step. After continuing for 1600 time steps, the two 
solutions agreed at every point to three significant figures. 

Although no difficulties are encountered in starting the scheme in 
the linear case, we found a weak nonlinear instability while solving the 

Navier-Stokes equations. Its onset could be detected by the values of u 11 

H t & xi 

and (u ) slowly diverging. The problem was cured by setting (u ) = 

u 11 at the end of the first few time steps. The calculation of u 11 (which 
t t 

is not needed elsewhere) almost doubles the computing time for each time 
step where it is needed In a test calculation on a 16 x 16 x 16 mesh, 
we found it was sufficient to make this correction after the first four 
time steps after which no further evidence of instability was seen for the 
next 130 time steps. In the mam 64 x 64 x 64 calculation, we alter- 
nately turned the correction on for four time steps and off for four time 
steps throughout the problem. An examination of the skewness of the veloc- 
ity field as a function of time (Section 5.1) leads us to believe that the 
correction could have been turned off permanently after eight applications. 

2.6 Space Differencing 

The time- integration scheme has been presented m terms of the ordi- 
nary differential equation 
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It is equally applicable to the system of incompressible Navier-Stokes 
equations 


3t 3x ^ U i U j^ 
1 


3p . „2 

- + VV U . 

3x^ i 


if the spatial derivatives are calculated at least as accurately as the 
time derivatives. The method for calculating Su^/St and maintaining 
zero numerical divergence will be discussed in Chapter III. The spatial 
derivatives were approximated by fourth-order accurate spatial differenc- 
ing. This means that the equations we are actually solving are 

+ —— (u u ) = - + W^u + 0(<5t 3 ) + O(Sx^) (2.26 

3 t 3x l 3 3x l 

1 i 


— = 0(<St 3 ) + 0(6x 4 ) 


To properly compare error terms, we must include the velocity c to be 
dimensionally consistent. When we compare terms of 0(6t n ) to terms of 
0(5x n ), we really need to look at terms like (c 6 t) n and 6 x n . The 
choice of the appropriate value of c to use m the nonlinear case is 
somewhat unclear in a turbulent flow simulation, but the r.m.s. velocity 
is probably a reasonable guess. However, the stability condition gives 
essentially an upper limit on cdt/dx, the Courant number, where now the 
maximum value of u must be used for c to assure safety. The result is 
that <$t is normally chosen so that the time-diff erencmg error is m 
fact somewhat smaller than the space-diff erencmg error. 

Now we consider conservation of momentum and energy. For simplicity 
we have chosen to use centered, fourth-order spatial differencing in space 
with all quantities cell-centered. The scheme which was used m the mam 
calculation was 


1 2 
-y(Duu + uDu)-Dp + vD 1 u 
2 j i J 3 3 i i ki 
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where 


D = a fourth-order approximation to 9/9x , 


D u 
3 i 


^(j-2) + 8[u i (j+l) - u^Cj-l)] - u^j+2) 


126x 


(2.28) 


and 


D, 


a fourth-order approximation to V : 


\ U i = + u^Ci-l) + u^j+1) + u^j-1) 

+ u^(k+l) + u^Ck-l)] - u^(i+2) - u^Ci-2) 

- u^Cj+2) - u^Cj-2) - u^Ck+2) - u^Ck-2) 

- 90 u }/126x 2 

i 


(2.29) 


In both cases the derivatives are approximated at (i,j,k) and obvious 
indices are suppressed. Kwak (1975) has shown that the term - (D^u + 
u D u ) is conservative of both momentum and energy This means that 
no momentum or energy are introduced as a result of the spatial differ- 
encing, i.e.. 


;(D ±Uj + u j D^u i ) = 0 


E u (D u + u D u ) = 0 

1 1 3 3 3 1 


where the summation is over all grid points. The method can be called 
semi-conservative, because some error is introduced by the time-integration 
scheme. 
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Chapter III 


THE MAIN CALCULATION 


3.1 The Basic Equations 

The purpose of the main calculation, is to obtain, as accurately as 
possible, the solution of the equations of motion for homogeneous iso- 
tropic turbulence in an incompressible fluid. Physically, this flow is 
produced by passing a uniform stream of fluid through a mesh to produce 
the turbulence and then observing its decay as the turbulence proceeds 
downstream Special care is necessary to assure isotropy, but the ex- 
periment has been successfully carried out several times; the most re- 
cent such experiment is that of Comte-Bellot and Corrsin (1971) . An 
alternative to grid turbulence is box turbulence, m which the fluid m 
a box is stirred up and allowed to return to rest. To simulate grid tur- 
bulence we will use the Navier-Stokes equations (3.1. a), together with 
the continuity equation for an incompressible fluid (3.1.b)‘ 


3u 

i 

at 


+ 


a 


3x 

3 


(vy 


Y 5 - + VV 2 U. 
3x. i 

i 



(3.1. a) 


(3.1.b) 


where we use the summation convention. 

We shall attempt to simulate the grid turbulence experiment by se- 
lecting a cube of fluid and following its history as it passes downstream 
from the grid, i e., we are following it in a Lagrangian sense by invok- 
ing Taylor's hypothesis. In order to do this successfully, we must assure 
that the cube of fluid we select is large enough that all correlations 
are essentially zero at distances equal to the side of the box. From a 
practical point of view, this means that the box must be large compared 
to the integral scale of the turbulence. On the other hand, the box 
should also be small enough that, under the conditions of the experiment, 
no significant changes m important integral properties occur over a dis- 
tance equal to the size of the computational cube. If these conditions 
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are met, we may simulate the experiment by following the time history 
of the cube of fluid using periodic boundary conditions m all three 
spatial dimensions. 

Equation (3.1.b) can be replaced by an equation for the pressure 
p. We apply the divergence operator to Eq. (3.1. a) and note from Eq. 
(3.1.b) that 


~ 3u 

-3L_i = o 
at 8x 


and 


V 2 3u. 
i 

3x 


0 


We are then left with a Poisson equation for the pressure. 




3u 
i 

3x 

3 


(3.2) 


In general terms, the method of solution is to start with the veloc- 
ity field at time ndt, solve Eq. (3 2) for the pressure field at time 
ndt, then use Eq. (3 l.a) to find 3u^/3t at n<5t and advance the so- 
lution to time (n+l)6t using the method described in the previous 
chapter. Thus we insure continuity at time (n+l)6t by properly choosing 
the pressure at time n6t. 


3.2 Derivation of the Pressure Equation 

The Poisson equation (3.2) for the pressure is exact, no approxima- 
tions are involved in its derivation. It is, nonetheless, instructive to 
examine the origin of the pressure equation from a numerical viewpoint. 
The final expression we arrived at m Chapter II for the velocity field 
at time (n+l)6t was 


n+1 


n 

u + 

l 


r 1 ( *\n-l , 2 *.n 5 , * n+l~] 

SC L~ 12 ( "i>t *3 ( Vt + 12 (u !>t J <3 ' 3) 


Now we require that the numerical divergence of u be zero, i.e., 

n+1 1 

Du =0, where D denotes the fourth-order numerical difference 

ii i 

approximation to 3/ax^ defined by Eq. (2 28). 

Suppose that at this point m the calculation we have already evalu- 
ated (u ) n \ (u ) n , and (u ) n+ ^ but not p* . We have from the 
it it i 

Navier-Stokes equations 
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(3.4) 


f *\ n+1 



[' 


* * 
D.u, u 
3 i 3 


4- JU O *A* 

+ u D . u + vD u 
] ] lj k i 


* 


°l P 


where all missing time indices are assumed to be n+1. This leaves p 
as the only unknown. We now substitute Eq. (3.4) into Eq. (3.3), apply 
the numerical divergence operator to the result, and require that D^u^ 
be identically zero. This yields 


W > 


12 n 
-=- Du + 
5 ii 


D i [- i 


-1,8. *.n 
+ 5 <u ih 


1 4- -> Ju JU 

ir (D u u + u D u ) + vD! 

2 1 i 1 111 


> * 
'u 

C 1_ 


(3.5) 


If Eq. (3.5) is solved exactly, the numerical divergence of u will 

1 2 

be identically zero (within computer round-off error) . The operator D^ 

m Eq (3.5) is the fourth-order numerical difference approximation to 

A 

the Laplacian defined by Eq. (2.29). The operation D (D p ) implies 

a 1 1 2 

two sequential operations of D on p . In one dimension, is a 

five-point operator. 


D 2 f(k) = [-£ (k-2) + 16f (k-1) - 30f (k) + 16f(k+l) - f (k+2) ] /126x 2 


and D D i is a nine-point operator given by 

Di( D i f (i)) = [f (i~4) - 16f (i-3) + 64f(i-2) + 16f(i-l) - 130f(i) 

+ 16f(i+l) + 64f (i+2) - 16f (i+3) + f (i+4) ] /1446x 2 

2 2 

Note that this is not the simplest fourth-order approximation to 3 / 3x^, 
but one must use this operator (and none other) to insure that the contin- 
uity equation is satisfied exactly. 

Depending on how the initial conditions are set up, the first two 
time steps may or may not have velocity fields whose divergences are 
identically zero. In our case the divergences were small but not zero. 

If we retain all of the terms m Eq. (3.5), including the non-zero diver- 
gences, the next two time steps will have identically zero divergences. 

As we explained in Chapter II, for the first few time steps of a calcula- 
n ti 

tion we set (u ) = u fc at the end of the time step for stability reasons. 
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All we normally have at the end of a cycle is u 11 , not u^. To calcu- 
late u n we need the pressure field' p n , which is not normally calcu- 

t ^ £ n 

lated. The pressure field p is calculated by requiring that D^u ) 

= 0. After doing this twice, we have made \ D^u^, 1 , 

and D (u,)^ all identically zero. By referring to Eq. (3.3), we see 

that the only additional requirement needed to make Du = 0 is that 

* n+1 

D^Cu^) =0. To satisfy this requirement, we apply the divergence opera- 
tor to Eq. (3.4) and we obtain 


D (D p ) = 

i i r 


- D 


n * * * * 

tt (D u.u + u D u ) 
2 3 i 3 3 J 1 . 


(3.6) 


which is considerably easier to work with than Eq. (3.5). Since the term 

-(Duu +uDu) forms a part of u , we do the actual calculations 
3 1 3 ] J i t 

in the following sequence. 

* n+1 


1. Calculate (u ) 
Chapter II. 

1 


using standard leap frog, as described m 


2. Calculate 


a A A 0 * 

„ (D u u + uDu)+D.u 
2 J i 3 2 3 i ki 


— , n+1 


and store the re- 
*. n+1 


suit in the disk file, which will later contain (u ) 

: k 

Using the result from step 2, solve Eq. (3.6) for p (Note 

2 ^ 

that the inclusion of D.u does not affect this calculation, 

k l 


since D 


2 * 

D. u = 0 . 

i |_ k i_ 


) 


Evaluate -D p and add the result to the results of step 2, 
1 * n+1 

thus leaving (u ) on disk. 


5 . Calculate u' 


n+1 


from Eq. (3.3). We have now completed one time 


step. 


3.3 Solution of the Pressure Equation 

The Poisson equation for the pressure is solved with the help of dis- 
crete Fourier series. Any one dimensional set of N 'numbers which repre- 
sent the values of a function f at N evenly spaced grid points x = 

3 Ax, 3 = -N/2 ... N/2 - 1, can be uniquely represented by a discrete 
Fourier series, i.e.. 


23 



K = l~l 


| E F(k) 


The Fourier coefficients are given by 

_ N , 

J 2 “ 1 


FCk) = E f(j) exp^jkJ 


N ' / 

J= ~2 

In these equations, f(j) represents the value of the function f(x) at 
the point x = j<$x. Likewise, let p(j) represent the value of the un- 
known function p(x) at x = jdx. Then the solution to the discretized 
Poisson equation is 


DDp(x) = f(x) 

X X 


(3.7) 


where 


n - , V = f(3-2) - 8f (j-1) + 8f(j+l) - f(j+2) 

x 126x 


can be found by substituting discrete Fourier series for p(j) and f(j) 
into Eq. (3.7): 


D x D x I>00 ex p(~^p J* k ^ = E F ^ k) exp(^i jkj 


(3.8) 


POO = EP ( 3) exp^— 1 jk^ 


It is easy to show that 


D D exp 
xx 


) = |E p 00 exp^-l^jkJ 
(“IP 3 k ) = ' 8 ( k > ex p( Z lf i J k ) 


(3.9) 
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where 


g(k) = - 130 + 32cos 


(¥) + - »~(¥) + 2 “ S (¥) 


Applying (3.9) to (3.8), we get from the linear independence of the com- 
plex exponentials 


g(k)P(k) exp(-|p jk 


Now we multiply Eq. (3.10) by exp 

k = | - 1 


)- 
(¥*■)• 


F(k) exp 




(3.10) 


and, noting that for k ^ k' 


k - - 


2 exp J(k-k’) I = 0 

N 1 — —1 


which leaves us with 


P(k) 


_ F(k) 
g(k) 


(3.11) 


With the above m mind, we see that the equation 


d x d x p(j) = f(l) 


can be solved by the following three steps. 


1. Transform f(j), i.e., compute F(k)=^f(j) exp^p jkj 


2. Calculate P(k) = P~- 

g(k) 


3. Invert P(k), i.e. , compute p(j) = |-X)P(k) exp JJ k j 

The extension of the method to three dimensions is straightforward. 
The solution to the equation 


is obtained as follows: 


1. Transform f (j » j 2 ^ 3 ), i.e , compute 


F( k l» k 2 ’ k 3^ = S 23 23 f(l 1 >J 2 ‘ J 3) exp P^ (3i k i + i2 k 2 + J 3 k 3M 

7 11 < — l 


3^ 3 2 ^3 
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2 . 


Calculate pflc^k^kg) = F(k 1 ,k 2 ,k 3 )/g(k 1 ,k 2 ,k 3 ) 

3. Invert Pfk^jk^k^), i.e. , compute 

pCj-plo.jJ = £ £ p ( k i > k o> k Q) 

- N k x k“ 2 k 3 

The use of the fast Fourier transform (FFT) algorithm (Cochran, 

1967), which requires CN log 2 N operations to perform the one-dimensional 
Fourier transform of N data points, makes the above method of solution 
practical. 

3.4 Modifications to Reduce the Running Time 

ft 

(a) Solving a Two-Dimensional Poisson Equation Using a One- 
Dimensional Fourier Transform 

The two-dimensional Fourier transform of data on an N x n grid is 
normally accomplished by performing N one-dimensional transforms m one 
direction followed by N one-dimensional transforms in the second direc- 
tion, for a total of 2N one-dimensional transforms each of length N. 

Since each one-dimensional transform of length N requires CN log„ N 

2 ^ 

operations, the method just described requires a total of 2CN 1 og 2 N 
operations. The constant C represents four multiplications and two 

additions. Suppose we could do the same thing with a single one-dimensional 

2 2 2 2 
transform on N points. This would require CN log 2 N = 2CN log 2 N 

operations, i.e. , exactly the same number as before. It turns out, how- 
ever, that the machine language fast Fourier transform routine used for 
our problem is twice as efficient (in computing time per point) in calcu- 
lating a 4096-pomt transform as it is in calculating a- 64-pomt transform. 
Thus, if we can solve the two-dimensional Poisson equation on a 64 x 64 
grid by using a single 4096-point transform, we can reduce the running 
time for that part of the problem by 50%^ 

In order to see how we might take advantage of this, consider the 
4x4 grids illustrated m Fig. 3 1. The points m parentheses represent 
the virtual points which are used to provide the boundary conditions. The 
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Fig. 3.1. a. Normal periodic boundary condition 
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Fig. 3.1.b. Modified periodic boundary condition 
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periodic boundary conditions which are normally used are illustrated m 
Fig 3.1. a. The required virtual data on each line are taken from the 
opposite end of the same line. A modified or staggered periodic arrange- 
ment is shown m Fig. 3.1.b, here the virtual data on the horizontal 
lines are taken from the opposite end of the succeeding or prior line. 

The only difference m the two cases is that the modified conditions of 

Fig. 3.1.b have the left and right boundaries offset vertically by one 

2 

cell. The arrangement of Fig. 3.1.b allows use of the single N point 
transform and is therefore desirable. This raises the question of what 
we should require of the boundary conditions for the box turbulence prob- 
lem. The first requirement is that the virtual data must have no correla- 
tion to the adjacent points. This means that the correlation of the veloc- 
ity field across the box must be negligible (assuming the virtual points 
are taken from the opposite end of the box) . This requirement is met since, 
as noted earlier, the box- turbulence problem must have a velocity field 
which is uncorrelated half way across the box if it is to make physical 
sense. Secondly, the row of virtual points must represent turbulence. 

This requirement can be met by equating the data at the virtual points to 
those of some other row of points m the problem. As long as the order of 
the points is retained (i.e., the statistics m the vertical direction are 
unchanged), they may be shifted vertically with no effect. Either of the 
boundary conditions illustrated m Fig. 3.1 satisfies the above require- 
ments The advantage of the boundary condition illustrated in Fig. 3.1.b 
will soon be apparent. 

Suppose we combine the four rows m Fig. 3.1 into one row of 16 points 
numbered j = 1,2,3, ... ,15,16. For purposes of the Fourier transform, it 
is more convenient to use an index m - j - 9 so that: 

m=7 /« \ 

F(k) = V f(m) exp(-^ mk) ( N = 4) 
m = - 8 > 

f(m) = ^F(k) exp (~y L ' mk) ( N = 4) 

N k V N ' 

Now we let j = j^ + Nj^. ^ en f (j -j+1 ,J 2 ^ = for a H points on the 

grid. Referring back to our derivation of g(k) in the preceding section, 
we see that if x = j^ 
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where 


D x D x f (m) = F(k) exp/' - 2 - ^ 1 mkj 

N v ' N > 


gl (k) = - 130 + 32cos^^j + 128cos^pj - 32cos^~j + 2cos^-~j 

f(j 1 >J 2 + 1 ) can be written f(j+N), with N = 4 m this case, and it 
immediately follows that if y = 32 ^y> 

D D f(m) = -^-^g„(k)F(k) exp(- mk) 
y y N k ' H ' 


where 


g^Ck) = - 130 + 32cos 


(“) 


+ 128cos 


32cos 


/ 6irk\ , „ / 8irk \ 

\“r ) + 2oos \tt j 


Thus, wherever N occurs m g^Ck), it becomes N m g 2 (k), since a 
difference of one grid point m the x direction is the same as a differ- 
ence of one unit m the j, but a difference of one grid point in the y 
direction is the same as a difference of N units m j. Hence, 

AT £ 


N 

k = T-l 


N 

k = 'T _1 


(DD+DD) E P(k)exp(-^£mk) = E„ gl (k)+g,(k) 
yy , N 2 ' N / , N 2 •- 1 Z -i 


k = - ; 


P(k) exp 


(1M 


and the two-dimensional Poisson equation. 


(D x D x +1D y D y) P(X ’ y) = f(X ’ y> 


is equivalent to 


N 

k = T-i. 


( D x D x + D vV ^ P< k ) ex p(~^T L mk ) = S „ f (k) exp/-™- mk) 
x X y y , _ n2 ' , N 2 \ N 2 / 


k= - 


and can be solved by the following steps: 


k =_ 

2 


1 F(k) =Ef( m ) exp/^y 1 mkj , where m = j, +Nj„- (N 2 /2 + 1) 
m \ N / 
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2. P<10 = F(k)/t 8 l (k)+g 2 (k)] 

3. p(m) = -^-^P(k) expf— 5^ mk 

N k ' N 

Exactly the same methods could be applied to reducing a three- 
dimensional transform to a one -dimensional transform, but the CDC 7600 
large-core memory is not large enough to hold all of the necessary data, 
and the full advantage of this method is unattainable. The three- 
dimensional transform is therefore done by a series of modified two- 
dimensional transforms on each plane of data followed by a one-dimensional 
transform m the third direction. This reduces the running time of the 
Poisson solver by one-third (a 50% savings on two thirds of the trans- 
forms) . Since the Poisson solver takes roughly half of the total running 
time, this makes the overall saving about 16% of the total. 

(b) Savings from Simplified Indexing 

The modified boundary conditions described m the preceding section 
have the additional advantage of allowing us to write all of our differ- 
encing equations m terms of one-dimensional arrays. To illustrate this, 
we show how a two-dimensional case can be reduced to a one-dimensional 
problem. Suppose we have a 64 x 64 array, u(64,64), and we wish to 
calculate its Laplacian to second-order at each point. The easiest way 
to program this would be to define a new array v(64,66) with v(i,l) = 
u(i,64), v ( i , j ) = u(i,j-l) for j = 2, ..., 65, and v(i, 66 ) = u(i,l); 
each of these relations holding for l = 1, 64. This takes care of 

periodicity in the 3 direction, and we could then write our FORTRAN 
program as 

do 10 i=l, 64 
ip l=i+l 

if (ipl. eq. 65) ipl=l 
iml=i-l 

if (iml. eq. 0) iml=64 
do 10 3=2,65 

10 ulap(i, 3 -l)=v(ipl, 3 )+v(iml, 3 )+v(i, 3 +l)+v(i, 3 '-l)- 4 .*v(i, 3 ) 
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There are two difficulties wxth the above codxng Firstly, the 
"if" statements used to calculate xpl and xml slow the executxon of 
the loop- Wxth the modxfxed perxodxc boundary condxtxons, we can sxmply 
wrxte x— 1 for xml and x+1 for xpl and we wxll be usxng the cor- 
rect poxnts even at the ends of the rows, thus elxmxnatxng the "xf" 
statements. Thxs reduces the runnxng txme of a typxcal loop by 10%. 
Secondly, txme xs requxred by the computer to fxnd the address of the 
varxable u(i,j). To fxnd thxs address xt must calculate m = i+64*j . 
Sxmxlarly, xn the case of a three dxmensxonal array, to fxnd the address 
of u(i,j,k) xt must calculate m=i+64*j+4096*k. Wxth our modxfxed 
boundary condxtxons, xt xs possxble to code the above loop as follows. 

do 10 m=l, 4096 

10 ulap (m)=v(m+129 )+v(m+127)+v (m+192)+v (m+64) -4 . *v (m+128) 

Use of sxngle subscrxptxng resulted xn an additional 30% savxngs xn the 
typxcal loop. Thus, the use of staggered perxodxc boundary condxtxons 
allows a reduction xn the runnxng txme of all the dxfferencxng calcula- 
tions (x.e., vxrtually all the calculatxons other than the fast Fourier 
transforms) by a net 40%. In fact, xt was thxs savxngs whxch prompted 
the xnvestxgatxon of the modxfxed Poxsson solver 

The combxnation of the modxfied Poisson solver, the sxngly subscrip- 
ted arrays, and the third-order txme scheme which allowed an increased 
time step reduced the total running txme of the main calculation by a 
factor of six. Whereas we initially anticipated the problem would use 
nine hours of computer time, the final version used only 90 minutes. 
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Chapter IV 


TURBULENCE MODELING 


4. 1 The Equations of Motion 


For an incompressible fluid, the Navier-Stokes equations (4.1) 
together with the continuity equation (4.2) describe the motion of a 
turbulent flow. These equations are 



U i U J 


_ 

9x, 


+ vV^u, 



(4.1) 

(4.2) 


The term (8/9x )u u in Eq. (4.1) accounts for the change in u. 

□ 1 J ^ £ 

at a point in space due to the advection of momentum. The term vV u_^ 

accounts for the change in u^ at a point in space due to viscous 

forces . 

Consider an eddy of size L whose typical velocity is U. For 

2 

such an eddy, the advective term is of order U / L and the dissipa- 

2 

tive term is of order vU/L . The ratio of the advective term to the 
viscous term is of order UL/v,, which is the (non-dimensional) Reynolds 
number Re. An eddy with Re << 1 is dominated by viscous dissipation 
and will rapidly die out. An eddy with Re » 1 is dominated by ad- 
vection and will remain in the flow for a relatively long period of time 
before it will die out. Hence, Re = 1 gives an estimate of the small- 
est eddy one would expect to find in the flow. We assign to the small- 
est eddy the size q and the velocity u^, and since Re = 1 for this 
eddy 


nu^ = v (4.3) 

In order to get another relationship between the velocity and 
length scales of the smallest eddy, we multiply Eq. (4.1) by u^ which 
gives the equation for the kinetic energy 
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(4.4) 


at 


U i U i 



u 


3 



- u 


_a 

3x 


+ u.vV u. 
i l 


Integrating Eq. (4.4) over a volume V, applying the incompressibil- 
ity requirement (4.2), and ignoring contributions from the boundary 
yields 



dV' 


v 



u v 2 u.dv' 

l i 


(4.5) 


The right-hand side of Eq. (4.5), when integrated by parts, is seen to 
be negative definite and thus represents the energy-dissipation rate. 
Letting e be this energy dissipation per unit volume and noting that 
the dissipation occurs mainly in the eddies of size q , we have 

2 

u v 

e - — y- (4.6) 

n 


Now, since the energy dissipated by the small eddies comes from the 
-large eddies, the dissipation rate e is really determined by the large 
eddies and can be regarded as given. Then the only unknowns in Eqs. 
(4.3) and (4.6) are the turbulent microscales q and u^. Solving for 
these we get the Kolmogoroff (1941) expressions for the turbulent mi- 
croscales . 


q = 




(4.7) 


A more detailed discussion of the Kolmogoroff microscales can be found 
in Tennelces and Lumley (1972). 

In a typical problem involving turbulence, the length scale of the 

largest eddies which we want to simulate is several orders of magnitude 

larger than q. In fact, L/q, where L is the largest scale m the 

3/4 3/4 

problem, is of order Re . Hence we would need Re grid points 

L L | 

m each direction, but the largest number of grid points m each dimen- 
sion which one could squeeze into present-day computers for a three- 
dimensional calculation is of the order of 100. On the other hand, 

typical Reynolds numbers of engineering and scientific interest are in 
4 8 

the range 10 -10 . Consequently, the grid point separation is, m 
general, orders of magnitude larger than p. Hence the eddies of size 
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"n, in which the energy dissipation is occurring, cannot be calculated 
directly. We attempt to resolve this problem by defining a new velocity 
field where the overbar denotes some sort of averaging process. 

It could denote a temporal average, a space average, or an ensemble 
average over many realizations. We then define u| by u^ - u^ + u^. 

Leonard (1973) has suggested that the appropriate averaging proc- 
ess for large-eddy simulation should be a local spatial average of the 
form 

u^Cx) = — J G(x-x' )u i (x’ )dV' (4.8) 


where V is a volume surrounding the point x over which u^ is to 
be averaged and G is a weighting function as yet unspecified. This 
process may be called filtering, as its effect is to remove the small- 

Uj the filtered 

1 “i 

ox large-scale field and u^ the subgrid scale field. 

The simplest averaging operation is to let G = 1 and V be the 
cubic volume with sides of length A whose center is at x. Then 


scale fluctuations from u^ m forming u_^. We call 


. Aa , Aa , Aa 
X 1 + T X 2 + T X 3 + T 




'73 / / / <‘ 1 (*r x i’ x 2' x 2’ x 3‘ x 3 )dx i dx 2 dx 3 


a _Aa. Aa Aa 
X l” 2 X 2~ 2 X 3 2 


(4.9) 


Unless otherwise noted we will take Eq. (4.9) to be the definition of 
u^ throughout the remainder of this dissertation. 

We now take Eqs. (4.1) and (4.2) and obtain their filtered counter- 
parts by multiplying each by the weighting function G = 1 and integ- 
rating over the cubic volume V to obtain 


at 


+ u.u 

3x i i 
3 J 


3p , _2- 

- i + vV u. 
3x, i 



(4.10) 


(4.11) 


We now make the substitution u 

l 

term and obtain 


u. + u' in the nonlinear advective 
i l 
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(4.12) 


where 



3 

3x 

3 


u.u 
l J 



2— 
W u. 



n . . = u u' + u’u. + u'u' (4.13) 

i 3 i 3 i 3 i 3 

is the subgrid scale counterpart of the Reynolds stress . 

We stress the fact that Eq. (4.12) is exact. We have defined new 
variables, but so far we have made no approximations. We also point out 
that u and u' are continuous variables defined at all points in 

ii 

space and tune, and they are m no way tied to the finite grid of points, 
which will be introduced later. 

4.2 Approximations to Solve the Filtered Navier-Stokes Equations 

In order to solve Eqs. (4.12) it is now necessary to make some ap- 
proximations. The testing of these approximations is the purpose of 
our main numerical simulation. The three most common approximations 
used are 


u u « u u (4.14) 

1 3 i j 


u.u' + u'u - 0 (4.15) 

i 3 1 3 

u'u' = f(u ,u ) (4.16) 

3-3 i J 

One of our major purposes m this work is to investigate these approxi- 
mations, test their validity, and suggest improvements. Numerical 
tests are given m Chapter V. We give a discussion of each of these 
approximations below. 


4.3 The Approximation (4,14) u^u = u^u 

Leonard has shown that Eq. (4.14) is probably a poor approximation 
m a turbulent flow. Consider a function f (x) defined m some region 
of space. If f is fairly smooth, we can approximate f(x) locally 
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by its Taylor series expansion about the point x , which is taken to 

~o 

be the center of the filter volume, V, over which f will be inte- 
grated to obtain f(x Q ). Substituting the expansion of f into Eq. 
(4.9), we obtain. 


f( ?o> ■ f (?o> + 2f 


4- 0(A 4 ) 
x a 

~o 


(4.17) 


where A a is the length of one side of the cubic volume V. Letting 
f = u x u j> we ha ve immediately the Leonard approximation 


L ij U i U j u i u j 24 8x k 3x k ^ U i U j^ 


(4.18) 


The last term in Eq. (4.18) will henceforth be referred to as the 
Leonard term. Ther e can be little doubt that Eq. (4.18) is a better 
approximation than u^u_^ “ u i u j • We now as k: What is the magnitude of 

the Leonard term in relation to the other terms in Eq. (4.12)’ 

We can get some idea of its size by considering the simple, one- 
dimensional Fourier wave u = exp(ikx) and the lindar filter of this 
wave defined by 


u(x) 



(4.19) 


from which it is easy to show that for u = exp (lkx) : 


u(x)u(x) 


sin 2 | 

'kA \ 
/) 

sin(kA ) 
a 

/kA ' 

\ 2 

kA 


) 

a 


u 2 (x) 


sin kA 

a 


kA 

a 


u u (4.20) 


We can now quickly test the accuracy of approximation (4.18). If we 
2 2 

approximate 9 /3x by a second-order space-centered finite difference 
on a grid with spacing A and apply it to a Fourier wave, we have 

o 


g2 [cos(2kA ) - 1] 

— n (uu) = 5 s u u (4.21) 

9x A 

g 
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We are interested in the ratio a = u u/u u. The exact value is 

seen from Eq. (4.20) to be (kA ) -1 sin kA . In Table 4.1 we have 

a a 

given a for various values of kA . The first column shows the exact 

a 

result. The second column gives the approximation (4.14) for which « 
is always unity. The third column is the result obtained from Eq. 

(4.18) if exact (Fourier) differentiation is used, while the fourth col- 
umn is the result obtained by using second-order finite differences, 

Eq. (4.21), with a grid spacing equal to A /2 , a value which will 

a 

later be shown to be appropriate. 

Table 4.1 
The Leonard Term 


kA 

a 

Exact 
Eq. (4.20) 

No Leonard Term 
(Eq. (4.14) 

Fourier 
Eq. (4.18) 

Second-Order 
Eq. (4.21) 

Gaussian 

0 

1 . 

1 . 

1 . 

1. 

1. 

m/4 

0.9003 

1 . 

0.8972 

0.9024 

.9023 

m/2 

0. 6366 

1 . 

0.5888 

0.'6667 

.6628 

7T 

0.0000 

1 . 

-0.6449 

0.3333 

.1930 

3m/2 

-0.2122 

1 . 

-2.7011 

0.6667 

.0247 

2tt 

0.0000 

1 . 

-5.5791 

1.0000 

.0014 


We see from the table that for kA = ir/4, a wave whose wavelength 

cl 

is eight times the averaging volume or sixteen times the grid spacing, 
the effect of the Leonard term is approximately ten percent and the ef- 
fect increases at shorter wavelengths. Waves with kA > it are poorly 

cL 

treated by any of the approximations. However, the filter removes most 
of these waves (as it must to avoid the numerical problem of aliasing) , 
so the problem is not as severe as it might seem. The oscillatory beha- 
vior of the filter at high wave number has caused some workers to re- 
place it with a Gaussian filter. Its values are shown in the last column 
and we see that Eq. (4.21) does an excellent job of matching it. 

We emphasize that the error being discussed here is not related to 
"subgrid scale turbulence," but arises from incorrect handling of the 
interaction between waves which are supported by the grid. 
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4.4 The Approximation (4.15) u_^f = 0 
(a) A Model for the Cross Term 

We again consider the one-dimensional wave u(x) = exp(ikx) and 
the filter defined by Eq. (4.19). It is easy to show that 


uu 


sin (~2 £ ) 


kA 


1 - 



sm kA 2 
a u 

kA 


(4.22) 


We will now develop a model for the term u 1 u j • From the definition 
u - u . Using the expression (4.17) for u , 

'•L, 


of u' we have u' = 
3 3 


u . . - VVj + 0(4’) 


(4.23) 


which implies 


u u 
1 3 


a — 2 4 

zl V U 3 + 0(A a } 


(4.24) 


Now we substitute u = 


u - 


u u' 
i 3 


u into 
3 

A 2 
a 


Eq (4.24) 


24 


( u V 2 u + u V 2 u* ) 
\ i 3 1 3/ 


(4.25) 


The use of Eq. (4.17) to obtain Eq. (4.23) assumes that is "fairly 


smooth". This will be true if u is reasonably close to u , i.e., if 

3 3 

i 

3 

1 which is nearly resolvable on the grid and not 


u; does not fluctuate too rapidly. This implies that we are only model 


mg that portion of u 

that portion which is entirely subgrid scale. 

2 

. Since V u’ fluctuates rapidly throughout the averaging volume and 

J 2— 
has a mean value of approximately zero, whereas V u is rela ti vely con — 

d _ 2 — 2 — 

stant throughout the averaging volume, we expect that uVu' <<uVu, 

-t 3 r j 

and we can neglect the last term in Eq. (4 25). The lowest-order approxi- 
— 2 — — 2 — 

matron to u V u is just u 7 u , so the lowest-order approximation to 
r 3 1 3 

Eq. (4.25) is 

4 — 4 a- 2- 

= urn’ = — "XT' u V u (4.26) 


C. 

13 


24 
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Eq. (4.26) is our model for the cross term. Clearly there is no physics 
built into this model. 

In Table 4 2 we compare the values of the cross term m the same 
manner as we did for the Leonard term m Table 4.1. The values given m 
the table are u u'/u u. In comparing the magnitudes of the Leonard and 
cross terms, we should recall, first, that for the values m Table 4.1 
it is the difference from unity that is important, and second, that the 
cross term will appear with a coefficient of two in the equation. Thus 
we see that for kA^ = ir/4 the cross term has approximately half the im- 
portance of the Leonard term. As a function of k, the cross term 
increases in size and then decreases. The approximation (4.26) for the 
cross term is not as accurate as the corresponding approximation for the 
Leonard term. These conclusions will be borne out by the results presen- 
ted m Chapter V 


Table 4.2 
The Cross Term 


Exact Fourier Second Order 


tv hi 

a 

Eq. (4.22) 

Eq. (4.26) 

(Eq. (4 26) 

Gaussian 

0 

0. 

0. 

0. 

0. 

tt/4 

0.0236 

0.0257 

0.0254 

.0235 

ir/2 

0.0705 

0.1028 

0 0976 

.0718 

7T 

0. 

0 4112 

0.3333 

.0982 

3u/2 

-0.4949 

0.9253 

0.5690 

.0376 

2tt 

0. 

1.6449 

0 6667 

.0058 


(b) Leonard and Cross Term Energy Dissipation 

Assuming that the model given by Eq. (4 18) is correct, the energy 
dissipation due to the Leonard term is given by 


e 

L 



V 2 (u u )dV ? 
i J 


(4.27) 
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Carrying through the differentiation 9/9x^ and noting that 9u^/9x^ = 0, 
we obtain 


A r _ - a /_ 3u \ 

"L * " 24 ) v u i 9x k 9x k ( U j 9 aT/ dV 


(4.28) 


Now we integrate by parts twice with respect to 9/3x k> giving 


A /* 9u _ 
a / — i .2— . 

= - T 7 / U T— A U C 

24 X 3 9x i 


(4.29) 


Leonard (1973) has shown that (4.29) can be approximated m the case of 
homogeneous isotropic turbulence by 


/a 3. 

35 2^/ 9u iV 

48 aNhij J 


(4.30) 


For the cross-term energy dissipation we can write 


A 2 a 2 

E c ■ I f jf "x + it f v \ it {v\) dv ' < 4 - 31 > 

The first integral in Eq. (4.31) is identically zero, since by carrying 


out the differentiation 9/9x we have 

3 


f u ~~ u V 2 u dV T = f u V 2 u ■— 

X 1 dx 3 1 3 X 1 3 3 x 3 


(4.32) 


where we have again used 3u^/9x^ = 0, and by integration by parts we 
have 


f u (u V 2 U ^dv* = - f u V 2 U dV' 

J v 1 9x -, \ 1 v X 1 J 3x j 


(4.33) 


Since the right-hand sides of Eqs. (4.32) and (4.33) are negatives of each 
other, 




(4.34) 


We take the second integral m Eq. (4,31) and integrate by parts to ob- 
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~ f U p(u V 2 u|dV' = - 2 I /* u V 2 u r^dV' 

24 -'v 1 9x u ' 3 V 24 ‘'v 3 1 3x j 


(4.35) 


Comparing Eq (4.35) to Eq. (4.29) we see that the cross-term energy dis- 
sipation and the Leonard term energy dissipation are the same. Further- 
more, since the skewness 

fa 1^/8 3 > 


<p 1 / Sx lj 2 >' 


is known to be negative, both terms remove energy from the flow. This 
will be verified by our numerical experiments. 

4 . 5 Models for the Subgrid Scale Reynolds Stress uha^ 

Having developed models, for the Leonard term and the cross term, we 
are left with the equations 


8 / N 3p , „ 2 — 8 , T 

— (u.u ) = -Tf+v?u -r-(L 


9t 3x^ i 3 


1 ax^ ij 13 13 


+ C + T ) (4 36) 


where 


(4.37) 



A 0 - 


A 

13 

= v (u iUj ) 

c 

13 

a 

" 24 

13 

‘ n ij - I \k 5 ij 


= u'u* 
1 3 


(u V 2 u + u V 2 u ) 

' i -1 1 -1 


- . 'kk 
P = P + — 


We have indicated that L and C_^ , while important, are really the 
result of interactions of the resolvable scale flow field with itself. 
Conversely, x is solely the result of the effects of the subgrid scale 
motions which cannot be resolved. Our only hope for modeling t is 
that the subgrid scale effects, averaged over ‘the filter volume, are some- 
how functions of the resolvable part of the flow, i.e. , the filtered vel- 
ocity field u . The most obvious condition that such a model must satisfy 
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is that, since it represents the energy transfer from the resolved large- 
eddy field to the small subgrid scale eddies which are dissipative, 

- / u — t dV' < 0 

J v 1 13 _ 

The simplest model which satisfies this requirement is the eddy viscosity 
model 


8 

3X J Tl3 


_ 2 — 
- v T V 


where the eddy viscosity v^, is an adjustable constant, 
we can model directly by 


T 


13 



Alternatively, 


(4.38) 


where 


could be a constant or a function of position 


appears m the momentum equation is (3/3x^)r_^ 


The term which 
which can be written 


3x 


v 


T\ 8x 


9u 


3x 


If is always positive, then this term can be shown to be dissipative. 

Since the time scale for the small-scale turbulence is much shorter than 
that for the resolvable scale, we expect that the small-scale eddies will 
adjust to the large-scale ones. It is therefore reasonable to suppose 
that the local subgrid scale Reynolds stress should be a function of the 
local level of resolvable scale turbulence. Pursuing this line of reason- 
ing, we let v m Eq. (4.38) depend on the local flow variables. We 

1 2-1 
require that it be positive and have units of (length) x (time) . The 

most popular such model, due to Smagormsky (1963), uses 


v 


T 


(cA a ) 2 



/ 3u 

3u \ 

II 

' 2 . 

1 

2 ' 

1 3x 

3x / 


V 3 

1 / 



(4.39) 


Until now the only way to verify a model such as Smagormsky ' s has 
been to observe its overall effect on the resolvable scale turbulence. 
With our numerical simulation we can directly compare t. with the 
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results of the model at each point m space. It will be seen later that, 

although the Smagormsky model of x is not as good as the models 

lj 

which have been developed for the Leonard term and the cross term, it is 
still reasonably accurate. We will also consider other possible models. 


4.6 Combining the Leonard Term and the Cross Term 


When the models for the Leonard term and the cross term are added 
together, we obtain 


L 

il 


+ C 


il 




(4 40) 


It is notable that m one dimension Eq. (4.40) reduces to 


L 

il 


+ C 


il 



(4.41) 


This expression is equivalent to the quadratic form of the artificial 
viscosity sometimes used m compressible flow calculations which was 
first proposed by Von Neumann (1950). The purpose of the artificial vis- 
cosity in compressible flow calculations is to smear a shock front over 
several cells. We note that this is precisely the effect of filtering 
the velocity field. If the field u has a step discontinuity, u will 
appear as a ramp of length 2A , which is exactly what the artificial 

cl 

viscosity attempts to do. A major difference is that we are proposing 
that this term be included everywhere m the calculation, whereas the 
traditional use of the artificial viscosity is only m regions of compres- 
sion. Furthermore, the present approach provides a more rational approach 
to the development of this model. Since m any flow calculation one can- 
not resolve detail smaller than one to two meshes , we believe that this 
term should always be included m a calculation. Of course, m flows 
with relatively small gradients, its effect will be small 
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Chapter V 


NUMERICAL RESULTS 

5 1 Results of the Main Calculation 

The purpose of the main calculation was to simulate a low Reynolds 
number turbulent flow field on a 64 x 64 x 64 mesh which represents, 
as accurately as possible, a realization of a true turbulent flow. The 
computed flow can then be considered as experimental data which can be 
used as input for the analysis of various schemes to model the effects 
of turbulence. In this section we will show that the computed field is 
in fact a good representation of real turbulence 

The experiment on which our simulation is based was reported by 
Comte-Bellot and Corrsm (1971) . The physical experiment was the measure- 
ment of the decay of grid-generated "isotropic" turbulence m a wind tun- 
nel A time history of the decay is obtained by employing the Taylor 
hypothesis This eliminates the mean flow field by assuming that the 
flow variables at two points m the wind tunnel separated by a distance 
L in the direction of the mean velocity U q are equivalent to the flow 
variables at two times separated by the time t = L/ U q of a flow with no 
mean velocity at the same point m space. 

The conditions chosen to be numerically simulated are given m Table 
5.1. U q is the mean flow velocity, 10 m/sec, and M is the size of the 
mesh which originally generated the turbulence, 2.54 cm. The initial con- 
ditions for the numerical simulation were set up to coincide with the data 

r 

at U^t/M = 240 The initial conditions were given the same total energy 
and energy spectrum as the experimental data, and a zero divergence, but 
were otherwise random The data are in the form of a one-dimensional 
energy spectrum, E^(k) > from which we computed the three-dimensional 
spectrum from the relationship (Batchelor (1953)). 

E(k) 2 k 2k |_k 2k E ll (k) 


(5 1) 
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Table 5 1 


Gross Propertxes of the Turbulent Flow 
U =10 m/sec, M = 2.54 cm 


■ 

J? 

Dxssxpatxon 

Rate 

Kolmogorov 

Mxcro-Scale 

Taylor 

Mxcro-Scale 

r a 

U T 
o 

11 

/ cm \ 

( cm ^) 

(cm) 

(cm) 


M 

Vsec 1 

v sec ' 

V 

240 

6 75 

145 

.069 

.845 

38 1 

385 

5.03 

48.5 

.091 

1 09 

36.6 
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where E(k) is the three-dimensional energy spectrum. The initial field 
does not represent true turbulence since it contains none of the local 
velocity correlations that exist m a physical field. It is these corre- 
lations which give rise to the subgrid scale Reynolds str_esses which we 
hope to model. We also note that the skewness, which is an indication of 
the presence of turbulence, is initially zero. The expectation is __ that 
as the equations of motion are integrated m time, a representation of a 
true turbulent flow will develop. 

Given the fixed number of mesh points in each direction, N, the 

physical size of the box of fluid must be determined. The box must be 

large enough that the velocity correlation at L/2 is negligible, and 

it must be small enough that the highest wavenumber k = N /L is large 

max 

enough to include essentially all of the energy dissipation spectrum. The 
size of the box was chosen to be 20 cm x 20 cm x 20 cm. Fig. 5.1 gives 
the experimental velocity correlation function R^(r^,0,0) where 


F 11 ( r 1 ,0,0) 


P u (r l ,0,0) 

P 1]L (0,0,0) 


(5.2) 


and 


eil (r l’°’ 0) = <u l (x l» x 2’ x 3 )u l (x l +r l’ x 2’ x 3 )> 

We see that a 20 cm length is sufficient to meet the condition that the 
correlation at L/2 be small. 

Figure 5.2 shows the three-dimensional energy spectrum E(k) and 

2 

the dissipation spectrum 2vk E(k) of the initial conditions. Since 

k = 10 cm \ we do indeed capture most of the energy and dissipation, 
max 

In Fig. 5.3 we show the three-dimensional energy spectrum E(k) , the 

dissipation spectrum D(k) , and the energy transfer spectrum T(k) of 

2 

the final numerical flow field D(k) is simply 2vk E(k). T(k) is 
calculated from 


~ E(k) = T(fc) + D(k) 


(5.3) 



was calculated using the numerical values of u^ 


and 
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Fig. 5.1. Velocity correlation function. 
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Fig. 5.3. Energy transfer and dissipation 


We now examine the results of the numerical simulation. We look 

first at the energy, as a function of time and its rate of change, the 

dissipation rate We are solving the exact Navier-Stokes equations with 

no model for turbulent energy dissipation. The only dissipation present 

2 -1 

m the equations is a result of the physical viscosity v = 0.14 cm sec 
If we have chosen a sufficiently fine grid the energy decay of the numeri- 
cal calculation will match the experimental data. As can be seen m Fig. 
5.4 this is the case With At = 0 0073 seconds, 50 time steps equal 
0 365 seconds, which is the elapsed time between U^t/M - 240 and 
U Q t/M = 385. The total energy m the numerical simulation at time step 
50 is 3 2% low. The dissipation rate, which is a more sensitive indica- 
tor, is 11 3% high. This is the result of too high a numerical transfer 

of energy from low to high wavenumbers. A shift to high wavenumbers will 

/ 2 

vk E(k)dk, more than the 

total energy, which is given by J*E(k)dk. 

So far, we have seen that our box is large enough to contain a sam- 
ple of fluid whose velocities are uncorrelated across the box and is small 
enough to calculate essentially all of the real dissipation The only 
question remaining is, "Has the flow field developed into a truly turbu- 
lent field 9 " The skewness of low Reynolds number wind tunnel turbulence 
has been shown experimentally to be approximately -0.4 (Batchelor (1953)). 
The skewness S is defined as 


S = 


< 


< 


' 3u 
3 

3x 


■au 

3 

3x 


> 


> 


(5.4) 


where < > indicates an ensemble average. The skewness of the numerical 
flow field, with the average taken to be the average over all of the grid 
points, is shown m Fig. 5 5. The skewness starts at zero, since the 
initial flow field does not represent true turbulence. The skewness plot 
indicates that after only 15 time steps we appear to have stabilized the 
skewness. The slight dips m the skewness which occur every eight time 
steps were mentioned in the discussion of the third-order time-diff erencmg 
scheme The third-order scheme has a weak instability which must be cor- 
rected for occasionally. In this run the correction was alternately 
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Dissipation Rate 




[T = .0073 N sec] 

Fig. 5.4. Dissipation rate and energy 
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Skewness 




turned on for four time steps then off for four time steps throughout 
the calculation. The periodic dips m the skewness coincide with the 
turning off of the correction. It appears that after about 30 time steps 
the correction was probably no longer needed. 

The tailing up of the skewness near the end of the problem is prob- 
ably due to the continuing accumulation of energy m the high wave num- 
bers. Looking at the skewness, we decided that time step n = 40 would 
probably be our best representaiton of true turbulence, and this time 
step was chosen for the analysis to be described in the following section. 

Another, much more convincing, argument that the flow is truly tur- 
bulent will be given m Section 5.4. 

5.2 Testing of Subgrid Scale Modeling 

Having completed the mam calculation, we now have a realization of 
a flow field which has the characteristics of physical turbulence The 
data, which we treat as if it were from a physical experiment, is given 
on a 64 x 64 x 64 mesh within a box which is 20 cm on a side. We now 
imagine placing a coarse 8x8x8 mesh over the physical space occu- 
pied by the original fine mesh, i.e. , each side of the coarse mesh is 
eight times a side of the original fine mesh The relation between the 
fine mesh and the coarse mesh is illustrated m Fig. 5.6. Within each 
cell of the coarse mesh we have the experimental value of the velocity ~ 
field u^ at 512 evenly spaced points. Now we need to know the value of 
the filtered velocity fluid at each point in the fine mesh. Recall 

that the filtered velocity field is a continuous function defined at all 
points m space and is independent of the definition of the coarse mesh. 
We use a simple box filter with sides of length 17A/8, where A is the 
mesh spacing of the coarse grid. In order to get an average at a point, 
we use the value at that point and an equal number of points on either 
side. This means, the number of points we sum over must be odd, hence 
17A/8 instead of 2A. The value of the filtered velocity component 
u^(i,j,k), where i,j,k are the coordinates of the point on the fine 
grid, is given by 

. i+8 j+8 k+8 

u (i, 3 ,k) = — x- EES u.(i‘ >j'»k') (5.5) 

* 17 J i ’= i-8 j ' =j -8 k'=k-8 
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Fxg. 5-6. Illustration of fine mesh inside coarse mesh 
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This equation is equivalent to a box filter with sides of length 17A/8 

where A is the mesh spacing of the coarse mesh. Calculations will also 

be made using an averaging volume with sides of length 9A/8. Having 

calculated u , we also have u! from its definition u = u = u' . 

i 1 111 

For illustration purposes we have randomly chosen a line of 64 points m 
the x^ direction and have plotted and u^(x^) for these 64 

points m Fig. 5 7. 

The remaining quantities m which we are interested are 


u u = X) X) 2 u (, u 

£ m 17 3 7* j ’ k' £ 1 


(5.6a) 




(5.6b) 




(5.6c) 


17 i’ j* k' 


We now restrict our attention to the quantities u^, u^u^ , U x U j ’ an( ^ 

u'u' at the centers of the 512 cells defined by the coarse mesh. The 
l j 

claim made for the turbulence models under investigation is that the vari- 
ables u u , If u' , and u’u' can be expressed as functionals of u We 
i J 1 J 3-3 1 

will now demonstrate the extent to which this is true for the case of low 
Reynolds number, isotropic, homogeneous turbulence of an incompressible 
fluid. 


5.3 The Energy Dissipation 

The equation for the kinetic energy of the filtered velocity field 
may be obtained by taking the scalar product of Eq. (4.3b) with u^ We 
obtain 


1 9 , , 

2 Ji ( \ u i> 


u (u u ) = 

i 3x i 3 


_ ^ ■ . ft - + vV u + 

i ax i i 


u i 9x ^ L ij + C ij + T ij^ 

(5.7) 


For our purposes the dissipation due to the real viscosity can be ignored 
since we find that it represents less than five percent of the total dis- 
sipation on the coarse grid. When Eq (5.2) is integrated over all space 
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u(i), u(i) cra/sec 
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Fig. 5.7. Sample of filtered and unf altered velocity fields 


the contributions from the pressure gradient and the nonlinear advective 
term vanish, leaving 


1 

2 


l 


* <\V dv ' 




C dv' 

(5 8) 
(5.8) 


Using our experimental values of L , C , and x for each point in 
the coarse mesh, we can calculate the dissipation rate due to each term. 



e 

i 


, 512 

512 2 
n=l 


u rT - L 

1J 


, 512 

312 E 
n=l 


u C 

i 3X 3 


512 


512 

E - 

n=l 


- 8 
U T 

i 8x ij 
J 


(5.9a) 


(5.9b) 


(5.9c) 


The results are plotted in Fig, 5.8 for e at four time steps, using 

the averaging volume of 9A/8 on a side. Some comments on Fig. 5.8 are 
m order. First, we note that at time step n = 1, i.e., the initial 

condition, the dissipation due to the subgrid scale Reynolds stresses 
is zero. This is what we expect, since the initial conditions do not 
represent true turbulence. As the flow develops and becomes more physi- 
cal, the dissipation term from the subgrid scale Reynolds stress rises 
rapidly before falling off as the energy of the flow decreases . This is the 
evidence referred to m Section 5.2 which indicates the flow has devel- 
oped a truly turbulent nature. The finite differences used to obtain 
Fig 5.8 were taken on the. coarse mesh to conserve computing time. It is 
more accurate to take differences of the filtered quantities on the fine 
mesh (recall that the filtered quantities are defined at all points m 
space), and this was done at the final time step to obtain the dissipa- 
tion from the cross and Leonard terms as well as the subgrid scale Rey- 
nolds stress. At time step 41 the dissipation rate due to the subgrrd 

2 3 

scale Reynolds stresses found in this way is 5.32 cm /sec , which is 
a bit less than that found by coarse mesh differencing, that from the 
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2 3 

cross term is 3.95 cm /sec , and that from the Leonard term is 2 82 
2 3 

cm /sec . We recall that the models predicted that the dissipation 
from the cross and Leonard terms should be equal. Though not equal, 
they are reasonably close. 

The dissipation for the large (17A/8) filter case was also calcu- 
lated at time step n - 41. The major difference from the previous case 
is a decrease m the dissipation due to the cross and Leonard terms. 

This is because increasing the size of the filter smooths the resultant 
field, causing a decrease m skewness and hence a decrease m cross- term 
dissipation. 


5 . 4 Correlations Between the Models and the Numerical Experiment 

We are now m the position of being able to make direct comparisons 
between the models for subgrid-scale turbulence terms and their numerical 
experimental values. We define the correlation coefficient C(M,X) be- 
tween the values of the model M and the experiment X as 


C(M,X) = 


< MX > 

<m 1 2 > 1/2 <x 2 > 1/2 


(5.10) 


where 


512 


<MX> = M(n)X(n) 


(5.11a) 


n=l 


< M > = 


1 512 9 

512 E « 2 <"> 

n=l 


(5 lib) 


< X > = 


1 512 9 

512 E X 
n=l 


(5.11c) 


If M and X are totally unrelated, then C(M,X) =0 If M is a con- 
stant multiple of X, i.e., if the model is exact, C(M,X) = 1. 

There are three levels at which comparisons can be made, and these 
correspond to how the terms appear m the equations For the moment we 
restrict ourselves to discussing the subgrid-scale Reynolds stress T 
The most direct comparison is at the tensor level, i.e , between the 
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experimental and modeled values of t However, the term which actu- 
ally occurs m the momentum equation is an acceleration vector 3 t /3x 

3 

We define the vector level of comparison to be that between the experimen- 
tal and modeled values of 3 t^/3x^. The scalar level of comparison re^ 

fers to the energy dissipation u (3x /3x ) produced at each cell m 

-*0 3 

the coarse mesh by the experimental t and the modeled . The pri- 
mary purpose of the subgrid-scale model is to remove kinetic energy at 
the correct portion of the flow, so the scalar level of comparison is 
important, and we will find it to be considerably better than the other 
two 


5.5 Tensor-Level Comparisons 
(a) The Leonard Term 
The Leonard term is defined as 


L 

^1 


ij 



kk 6 ij 



u u - u u 
1 3 1 3 


(5.12) - 


and the model we use is 


n 


1J 


a 




1_ 

3 


a kk 6 ij 


a 

ij 


A 

a 

24 


V 2 (u u ) 
i 3 


(15.13) 


Fourth-order differencing has been used in evaluating all of the models 
which we will be discussing. Fourth-order space differencing gives one 
to three percent better correlations than second-order differencing. The 
differencing was done on the coarse mesh, because this is the mesh which 
would be available in a real simulation. We note that we can do better 
m this case, since we have on the fine mesh as well as on the coarse 

mesh and can „get a better approximation of the actual derivatives. We 
compared the results obtained by differencing on the fine mesh to those 
obtained on the coarse mesh and found the differences to be minor (the 
correlations for the Leonard, cross, and Reynolds terms m the case of 
the small filter changed from 0.909, 0.790, and 0.277 to 0.934, 0.744, and 
0.297, respectively). 
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We find that the correlation between the model (5 13) and the ex- 
periment (5.12) is 0.935 for the large (17A/8) filter and 0.909 for 
the small (9A/8) filter. The ratios of the r.m.s. value of the model 
to the r.m.s. value of the experiment is 1.60 for the large filter and 
0.788 for the small filter; the reason why these values differ from 
each other and from unity are not understood. Since the model for the 
Leonard term involves the fewest approximations of the three terms we 
are considering, we expect it to be the best, which it is. Also, since 
the model is based on a Taylor series expansion of the filtered velocity 
field, we expect the smoother velocity field produced by the larger fil- 
ter to give better results, and it does. The correlations as functions 
of i,j are detailed m Table 5.2. Only small differences are observed 
throughout the table. 


(b ) The Cross Term 

The cross term is defined as 


C 

iJ 


il 



kk.^13 


n . = u.u' + u'u (5.14) 

il i 1 i 1 


and the model we use is 


M = a - -7 a. . 6 ; a 

ij 13 3 kk 13 13 


- 77 (u A +u V 2 u ) (5 15) 
24 1 3 3 1 


In this case we find that the correlation is better for the smaller fil- 
ter than for the larger filter. This is probably because the experimen- 
tal values are smaller for the large filter than for the small filter, 
due to the smoother flow field. The correlations of 0.685 and 0.790 are 
less than for the Leonard term, but the r.m.s, ratios of model to experi- 
ment are better, i.e., 1.23 and 0.96. Details are in Table 5.3. 


(c) The Subgrid-Scale Reynolds Stress 

The definition of the siibgrid-scale Reynolds stress is 

V. 


T 

IJ 




kk 6 ij 


13 


u’u’ 
1 J 


(5.16) 
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Table 5.2 


O' 

w 


Correlation Between Exact and Modeled Leonard Terms 


A 


A 


17A 

8 



9A 

8 



1 

2 

3 

N 

1 

2 

3 

1 

.92 

.94 

.94 

1 

90 

.91 

.93 

2 

.94 

.93 

.93 

2 

91 

.89 

.92 

3 

.94 

.93 

.94 

3 

.93 

.92 

.90 


Average = .935 


Average = .909 



> 


> 


1/2 

172 


< M 2 > 


< L 2 > 


1/2 

172 


1.60 


788 





Table 5.3 

Correlation Between Exact and Modeled Cross Term 


u> 


\i 1 2 3 

, 


.69 

.67 

.72 

67 

.64 

.70 

72 

.70 

.65 


\i 1 2 3 


00 

r-* 

76 

82 

.76 

78 

80 

82 

o 

00 

.78 


Average = . 685 


Average = .790 



> 


> 


1/2 

172 



> 1/2 
> 1/2 


1.23 


.96 





The four models we use are 


M 


a 


1 

3 


a 6 
kk 13 


a 


3-3 



where K is given by 


model 1 


model 2 


K = 


(GA ) 2 

a 


K 



- (CA ) 2 (to to ) 1/2 
a 11 



1/2 


(5.17) 


(5.18a) 


(5.18b) 


model 3 K = (CA^) T (u^.u^) (5 18c) 

model 4 K = C (5.18d) 

In Eq. (5.18b), represents the vorticity co^ = £ ijk^ li k^ X ^ All 

four models were found to be equally valid with the best correlations found 
to be 0.363 for the large filter and 0.303 for the small filter. In model 
3 the value for u^u^ is taken from the experimental data. Although these 
correlations are considerably below those for the cross term and the 
Leonard term, they are clearly significant. The constants m the models 
were obtained by matching the r.m.s. value of the exact quantity with that 
of the model prediction. Detailed comparisons for the Smagorinsky model 
are shown m Table 54, and a summary for all models is given m Table 5 5 


5 . 6 Vector-Level Comparison 

(a) The Leonard and Cross Terms 

In the previous section we compared the models directly with the cor- 
responding stress tensors. Here we make the comparison with the terms 
which actually enter the momentum equations, i.e. , 

_1_L . -i-c 

3 Xj L u ’ ^ So 
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Table 5.4 


Correlation of Subgrid-Scale Reynolds Stresses 
with Smagorinsky Model 



Average = .346 Average = 277 


C 


.269 


C 


.247 





Table 5 5 


Summary of Correlations between Exact Subgrid Scale 
Reynolds Stresses and Models 


Term 

Model 

m 


■31 

1 * 

R. 

ij 

Smagormsky (5.18a) 

.346 

.277 

270 

.247 


Vorticity (5.18b) 

.344 

.260 

.294 

.275 


T.K.E. (5.18c) 

.363 

.303 

.196 

.175 


Eddy viscosity (5.18d) 

.352 

.295 



IgH 






mm 

Smagormsky 

.425 

,346 

.240 

.264 

j 

Vorticity 

.408 

.327 

.220 

.247 


T.K.E. 

.434 

.362 

.138 

.155 


Eddy viscosity 

.426 

.356 



3R 







Smagormsky 

.710 

.580 

.186 

.171 

3 







Vorticity 

700 

.582 

.202 

.191 


T.K.E. 

.723 

.606 

.085 

.095 


Eddy viscosity 

716 

.605 
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The results for the large filter show that the correlations range from 
0.935 to 0.947 for the Leonard term and from 0.685 to 0.689 for the cross 
term. For all practical purposes these are the same as for the direct 
comparison. 


(b) The Subgrid-Scale Reynolds Stress 

In contrast to the case for the Leonard and cross terms, we find a 

significant increase m the correlations between 9 t /8x and its ex- 

13 1 

perimental value over the direct correlation between the stress tensors. 
The results shown in Table 5.5 show that all models are again equally 
good, but the correlation has typically increased from 0.35 to 0.42 for 
the large filter. Comparable increases are seen m the small filter re- 
sults. The reason for this increase is not understood. We note also 
that the model constants decrease, with one exception; again the reason 
is not understood. 


5 . 7 Scalar-Level Comparison 

(a) The Leonard and Cross Terms 

Here we make our comparisons based on the terms which enter the 
energy equation, i.e., u^ (SL^/Sx^) and u i (9C^/9x_^). Summaries 
for the three levels of comparison are given m Table 5 . We see a small 

decrease in the correlations from the vector to the scalar level for both 
the Leonard and the cross terms. The relatively large disagreement m the 
magnitude of the dissipation due to the Leonard term and model are not 
considered serious, since the dissipation due to the Leonard term is 
relatively small 

(b) The Subgrid-Scale Reynolds Stress 

We see a very sharp increase in the correlations for the subgrid- 
scale Reynolds stress at the scalar level. For example, at the vector 
level the Smagormsky model with A = 17A/8 had a correlation of 0.425, 
but at the scalar level it is 0.710. Part of the increase may be due to 
the fact that both the experimental and modeled terms have mean values 
which are significantly positive. Even so, when the mean values of both 
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are subtracted out, the correlation between the fluctuating components 
of the exact and model values is still 0.535. We also note a further 
decrease in the model constant. 

5.8 The Subgrid-Scale Eddy Coefficient 

The models (5.18) contain constants which are usually called the 

subgrid-scale eddy coefficient. The value of the constant has no effect 

on the correlation between model and experiment As mentioned above, we 

can, however, adjust the constant to match the r m.s. values of the model 

to experiment. The values of the constants found in this way are given 

m Table 5.5 and were mentioned earlier. The constants obtained decrease 

as we pass from the tensor level of comparison to the scalar level. 

Since the primary function of these models is to represent the transfer 

of energy from large to small scales, which acts like a dissipation m 

the large scales, we recommend that the values given for the scalar level 

of comparison be used. For the Smagormsky model, these values are m 

excellent agreement with theoretical and experimental values which range 

from around 0.13 to 0 21 (Deardorff (1971)). We note that when the Sma- 

2 

gormsky model is formulated using the term (CA ) , the value of C is 

cl 

nearly independent of A ; this would not be the case if the grid spac- 

3 . 

mg A were used. It is encouraging that we have obtained about the 
S 

same value for G as is obtained by theoretical arguments assuming an 
inertial sub-range and by numerical experiments at high Reynolds numbers, 
even though we are at low Reynolds number and have no discernible inertial 
range This leads us to speculate that C is relatively independent of 
the spectrum of turbulence, at least m the isotropic case The values we 
obtain are within ten percent of those found by Kwak et al. (1975) and 
Shaanan et al . (1975) by matching model calculations to experimental 
energy decay. Since a change in numerical method can result m a ten 
percent change m the constant, we can say that we have indeed predicted 
the model constant without reference to experiment 
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5.9 Comments on the Correlations 


A striking result that can be reached from looking at Table 5.5 is 

that all four of the proposed models are essentially equally valid. 

Since all of the models use a positive scalar times (3u /3x + 3u /3x ), 

i J J 

we checked to see how often the sign if t concided with the sign of 
(3u i /3x^ + Su^/Sx^) and found it to be only 68%. We also ran a calcula- 
tion with K being arbitrarily adjusted at each point in space so as to 
give the best possible correlation. At the tensor level of comparison, 
we achieved a correlation of 0.51, versus numbers around 0.35 for the 
models considered above. We conclude that no model of the form (5.17) can 
do significantly better than Smagonnsky's . This includes models which 
attempt to calculate the transport of turbulent kinetic energy and models 
which attempt to calculate both the turbulent kinetic energy and a length 
scale or the dissipation (so-called two-equation models) . This is par- 
tially verified by the results of model (5.18c), which show that even if 
one could calculate exactly the turbulent kinetic energy m each cell 
this would not give a significant improvement. 

We also include m Table 5.5, at the vector level, a modified Smago- 
nnsky model where, instead of 


3x 


3 L 


K 


' 3u 

i 

3x 


+ 


3u 

3x 


we used 


KV 2 u 


(5 19) 


The correlation decreased only slightly. The Smagor insky model has the 
disadvantage that its finite-difference form does not detect a wave with 
k = ir/A, i.e , a sawtooth, since its first derivative is always calcu- 
lated as zero. This can result m the failure to dissipate sufficient 
energy at high wave numbers. The modified Smagonnsky model does detect 
and dissipate these waves by using the finite-difference approximation 
to V . Model (5.19) has the disadvantage that one cannot rigorously 
prove that it is dissipative 
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5.10 Other Models, Which Were Discarded 

All of the models considered above are reasonably good We list 
here some of the more reasonable-looking models which were tried but dxs- 


carded. 

lations 

The following three t-ensor eddy viscosity models all had 
of less than 0.02 with the numerical experiment. 

corre- 



T = CA 2 D . D. 

13 a 1 k kj 

(5.20a) 


T 

13 

■ “ai 

(5.20b) 


T 

13 

- “al< D ik D kj + VW 

(5.20c) 


where D is the strain rate tensor, 
id 


and R 

ij 



is the rotation tensor, 



The next three models were proposed because of their similarity to the 

Leonard term, and all had correlations with t of less than 0.02. 

ij 


3 3 

u u 


ij 3^ 3x fc 1 j 


(5 21a) 


0 0 

T ij 3x 3x U k U k 

ij 


(5.21b) 


13 


3 3 


3 3 , 


l 3x 3x (u j U k } + 3x 3x ^i^J 
L 1 k J 3 k - 1 


(5 21c) 
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Chapter VI 


CONCLUSIONS AND RECOMMENDATIONS 


Most of the following conclusions are strictly valid only m the 
case of low Reynolds number homogeneous isotropic turbulence. For some 
of them, the range of validity extends beyond the range for which it has 
been proven; for others the validity of such extensions is unclear. 

Our first conclusions and recommendations are concerned with the 
simulation itself. 

1. With the present computer capacity it is possible to simulate 
homogeneous isotropic turbulence accurately m three dimensions. The 
limitation to a 64 * 64 * 64 grid restricts the Reynolds number based 
on Taylor microscale to R <40. 

A *■' 

2. The use of the third-order time method that we have developed 
allows a considerably greater time step to be used with very little sacri- 
fice m computational time or accuracy. We recommend the use of third- 

or fourth-order methods m future simulations, and some work should be 
done to find the optimum such method. 

3 The use of staggered periodic boundary conditions allows a con- 
siderable increase m computational efficiency at no cost whatsoever. 

4. The results of our simulation agree with the results of the cor- 
responding experiment m all significant statistical quantities, and we 
are confident that they may be used for model testing. 

5. With larger computers that will be available m a short time, it 
will be possible to use 256 * 256 x 256 grids and increase the Reynolds 
numbers considered by a factor of four. We believe that these computations 
are important and should be done 

The next set of conclusions and recommendations is concerned with the 


models used to represent the subgrid scale turbulence. 

6. The Leonard term is indeed of considerable importance m the pre- 
diction of turbulent flows and should be included m any simulation. The 


approximation to (u u - u u ) suggested by Leonard is fairly accurate, 

-t 1 J 

although some adjustment of the constant may be desirable. An alterna- 


tive is direct computation of the terms. 


71 



7 • The cross term, ( u 1 u ' + u ’ u )> which has been neglected by 
many previous authors, is also important, although less so than the 
Leonard term. We have suggested a model of this term which appears to 
do an adequate job of approximating it. 

8 Eddy viscosity models do only a fair job of matching the actual 
subgrid-scale Reynolds stresses, but tfqey do a better job m matching 
the acceleration produced by the stresses and they are rather good at 
predicting the dissipation or energy transfer to the small scales. 

9. All models of the eddy-viscosity class that we tested seem to be 
approximately equally good, and we are unable to choose among them on the 
basis of this study. The constant eddy-viscosity model is essentially 
what has been used by Orszag and co-workers, and our results partially 
explain their success. This point is probably closely related to the 
Reynolds-numb er independence of the large eddies. 

10. Further improvements in subgrid-scale modeling are not likely 
to result from attempts to find improved formulas for the eddy viscosity. 
We have shown that the best any eddy-viscosity model could do is a rela- 
tively small improvement on the Smagormsky model. Thus, turbulent ki- 
netic energy and two— equation models which have been popular methods for 
time-average modeling are not promising approaches to subgrid-scale 
modeling. 

11. We were unable to find improved models for the subgrid-scale 
Reynolds stresses, although a number of possibilities were tried. Further 
work m this direction could be fruitful. 

Other recommendations that we would like to make are: 

13. The effects of strain and shear on turbulence are very impor- 
tant, as they occur m essentially every flow of technological interest. 
The approach of this report ought to be extended to include those cases. 

14. The subgrid-scale Reynolds stresses play a role m large-eddy 
simulations similar to that which the usual Reynolds stresses play m 
time— or ensemble— average calculations. We therefore suspect that analo- 
gous models ought to be equally valid m the two cases. Further work is 
needed to substantiate this suspicion, but, should it prove to be the case, 
our work would have important consequences for turbulence modeling in gen- 
eral. 
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15. One can derive exact equations for the subgrid scale Reynolds 
stresses. Using the approach of the present report, we can evaluate all 
of the terms m this equation and thus determine their importance and 
examine methods of modeling them. 

16. The approach used m this report can be applied directly to the 
testing of tune— and ensemble-average models. If, we can compute flows m 
which the modeled effects are present, we could test the models in a man- 
ner similar to that used m the present work. 
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Appendix A 


Higher-Order Time Methods 

In Chapter 2 we developed a third-order time-advancement procedure. 
This method is a predictor-corrector method which requires only one eval- 
uation of the time derivative per time step. In this appendix we look 
at two families of related methods in a more general way. Since the equa- 
tions we deal with are parabolic with respect to the time variable, it 
is sufficient for initial study to consider only the ordinary differen- 
tial equation 


u = au (A.l) 

where the dot denotes time differentiation and a may be complex. 

The two most important questions relative to a numerical method are 
accuracy and stability. Accuracy is usually defined by assuming that if 
u(0) were known exactly then the computed value of u(A) , which we call 
u(A), is related to the exact value by 

u(A) - u(A) + const. A 11 "*"'*' u*' 11 "*^ (A. 2) 

where 1S th e (n+l)st derivative of u. A method with this 

fch 

property is called n order. Loosely defined, stability means that 
the computation does not blow up. One common definition is that when 
the method it applied to Eq. (A.l) with a having negative real part, 
it does not produce a growing solution. For most methods, stability de- 
pends on the size of the step chosen, i.e., it is conditionally stable . 

\ 

Implicit methods may be unconditionally stable, but they are difficult 
to apply to nonlinear problems. 

An approach to these questions is to look at the solutions of the 

ctt 

difference equations. Although Eq. (A.l) has only one solution e , 
the difference equations may have multiple solutions. One of these ap- 
proximates the solution of the differential equation with the desired 
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accuracy, the others are called parasitic solutions. All roots must 
have magnitude < 1 for stability. 

For our purposes we want a method with the following properties 

1. High accuracy — this allows a large time step with acceptable error. 

2. Stability — the method must be stable for time steps as big as nec- 
essary for the accuracy requirement. 

3. Few function evaluations — m partial differential equation solving, 
the "functions" are partial derivatives and are costly to evaluate. 

4. Minimum number of different values of variables required — in par- 
tial differential equation solving, the "variables" are large arrays 
which require considerable memory See Chapter 2 on this point. 

The popular Runge-Kutta method has the first two properties but not the 
last two. Essentially what we will do is accept poorer (but sufficient) 
stability m exchange for properties 3 and 4. 


Two Evaluation Methods 

The proposed methods are two-step (two previous values required) 
predictor-corrector type methods. The most general such method is 


u. = a.u + a 0 Au + a 0 u , + a.Au , , 

* In 2 n 3 n-1 4 n-1 


(A. 3) 


u n+l = 6 l U n + e 2 A % + 6 3 U n-l + B 4 Au n-l + B 5 U * + 3 6 Au * ‘ 


These can be combined to give 


u 


n+1 


. 2 

a„u + a. Au + a„A u + b n u , + 
On In 2 n 0 n=l 


• 2 

b.Au + b„A u .. 
1 n=l 2 n-1 


where 


(A. 4) 


a o = 

6 1 + a i B 5 

b o 

= e 3 + a 3 e 5 


II 

l—l 

P 2 + a 2 3 5 + « 1 3 6 

b l 

= 3 4 + 3 5 a 4 + 3 6 « 3 

(A. 5) 

a 2 

a 2 3 6 

b 2 

= “4 B 6 



Thus, although there are ten constants (0^,3^), only six combinations 
actually matter, four constants can be chosen arbitrarily to make the 
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method as simple as possible. By choosing the a,, properly, it is 
possible to obtain a fifth-order method, the method so obtained is highly 
unstable. So we will give up one order of accuracy to obtain stability. 
Applying the method (A. 3) to Eq. (A.l) and looking for solutions of the 
type = p , we find that p must be a root of the quadratic equation. 

2 2 9 

P - (a 2 (aA) + a^(aA) -+ a Q )p - (b 2 (aA) z + b-^aA) + b 0 ) = 0 (A. 6) 

For the method to be stable at all, i.e., for both roots of this equation 
to be smaller than unity as A -»• 0, we must have 


For minimal accuracy, i-e., that one root approach unity as A 0, we 
must have 


a 


0 



To obtain higher-order accuracy, we match the coefficients of the Taylor 

ctA 

senes of one of the roots to that of e and find: 


1st order 
2nd order 
3rd order 
4th order 
Solving, we have 


- a l + b o b l “ 1 

- 2a 2 = 2a^ + bO - 2b 2 = - 3 

6a 2 + 3 a;L - b Q = 7 

12a 2 + 4a^ - b^ = 15 


b 0 ~ 1 


a 2 


17 - jb 
12 


b 0 + 3 


b o + 7 
12 


By choosing values of b n within the allowed range, we can generate a 
family of methods. It turns out that b^ = 1 has the poorest stability 
properties, while b^ = -1 has the least accuracy 
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The stability bounds for selected values of b Q are given m Table 

A.l They were computed m the manner described in Chapter 2. We see 

that maximum stability is obtained at bg - .25, and the allowable time 

step is considerably below that of the fourth-order Runge Kutta method, 

for which [ ctA I =3 8. However, m our calculations the effective 
1 'max 

value of |aA] is approximately 0.2-0. 3 for accuracy reasons, so this 
stability limit causes no problem. 

Table A.l 


Stability Bounds 

b o 

|aA | 

' 'max 

-1. . 

0 

0 . 

51 

0.2 

.60 

0.3 

.60 

0.5 

.55 

1.0 

0 . 


We also solved Eq (A.l) using this method with b Q = 0, with two 
different sets of constants. Within roundoff error, both methods pro- 
duced identical results. It should be noted that, to minimize roundoff 
errors, one should choose sets of with the smallest values pos- 

sible. 


One-Evaluation Methods 

The methods described above require two evaluations of derivatives 
per step; i.e. , u A and u^ both need to be computed. To avoid this 
we would need either B^ = 0 or a 2 = = 8 2 = B^ = 0, either of which 

is incompatible with Eqs.(A.5). In order to obtain a single- evaluation 
method, we therefore use 


U *n+1 = “l U n + “2 Au *n + a 3 U n-l + a 4 Au *n-l 

U n+1 = 3 l u n + B 2 Au *n + B 3 u n-1 + B 4 Au *n-l + B 5 u *n+1 + B 6 Au *n+l 


(A 7) 
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Which requires the evaluation of u A only. We now assume solutions of 
the form 


hA 


l u A 

n 




P 

K> 


u 



and find that p satisfies a quartic equation 

p4 ~ ( n l + h 2 (aA)|p 3 - |ri 3 + n^(ctA)|p 2 + n 5 (aA)p + n 6 (aA) = 0 

7 7 (A. 8) 

where 


n l 

= *1 

+ 

B 5°l * 

^4 

= - a 2^1 + a 4 + a l®2 + a 3^6 ’ 

^2 

ii 

a 

w 

+ 

a l 3 6 5 

^4 

= a 2 e 3 + - ct.^ - , (A. 9) 

n 3 

3 3 

+ 

a 3 6 5 ’ 

^6 

= a 4 8 3 - a 3 g 4 . 


So, again, fifth-order is the maximum possible. The method so obtained 
is again unstable, so we must settle for fourth-order. We now find that 
for stability 


0 1 n l 1 2 * 


and accuracy requires that the other parameters be related to rt^ by. 


n 2 3 8 n l 


4 , 


24 "1 


n 3 1 n l 


1 . _1_ 

3 24 n l 


1 

3 


19 

24 n l 


The stability limits for these methods is shown m Table A. 2. Maximum 
stability is obtained with = 2, and the .allowable time step is only 
slightly smaller than that for the two-evaluation method 

In testing these methods, however, we found that these methods do 
not always produce the accuracy that one might expect. In particular, if 
an arbitrary version of this method is used rather poor results are ob- 
tained unless the method is started carefully. The problem can be cured 
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Table A. 2 


Stability Limits 


n. 

oA . 

1 

' ‘max 

0 . 

0 . 

1 . 

0.30 

1.5 

0.42 

2.0 

0.53 


by requiring that the predictor step be accurate. The predictor can be 
made third-order accurate, but only if = 0, which results m insta- 
bility. We therefore recommend that the method be used with second-order 
predictors. For these the method is uniquely defined by the predictor 
step. Two possibilities are (1) using the leap-frog method as a predic- 
tor (n 1 = 16/9) . 


u . , = u . + 2Au , , 
*n+l n-1 n+1 


(A. 10) 


U- " J ( 16u n " 7u *n+l ) + Yf (46 “*n " U "*n-1 + 13 “*n+l } 


n+1 


and (2) using Adams Bashforth as the predictor (n^ = 20/11): 


U *n+1 = u n + 2 (3u *n ~ u *n-l> 


U n+1 ~ U ( " 9u *n-l + 20u *n+l > + Y5 ( ~ 86 ' 1 *n +16 ' 1 *n-l + ^n+P 


(A. 11) 


Good results were obtained with both of these methods. However, we note 
that both methods contain some large coefficients in the corrector steps 
which is undesirable from the point of view of roundoff error propagation, 
and we have used the safer third-order method described in Chapter 2 in 
our calculation.. 
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Appendix B 


PROGRAMS AND FLOW CHART 


Flow Chart of Main Program 



Location m 
Listing 


Main. 1 -+• 1001 


Mam 1166-1541 


Main: 1543-2407 


8 ] 










Mam: 2*11-2555 
and subroutine 
DVDTMP 


Mam. 3533-4301 
and subroutines 
CALCPR and FFTFXY 


Mam: 4303-5053 
and subroutine 
DBLSQ 


Mam. 5251-6137 


Mam: 5251-6137 
and subroutines 
FFTBXY , DVDT 


Mam: 6141-10374 
and subroutine 
ADVNC 


Mam: 10376-11547 











( 1 ) 


For stability reasons we sometimes recalculate (u ) using the 
corrected third-order u 11 . 

(2) For the first three time steps only, the divergences of the previous 
time steps are not necessarily zero. 

System Subroutines Used m Mam and Start 

1. MEMREQ 

Call MEMREQ (NWORDS) request NWORDS of LCM space for this program. 

2. FORQTS 

Call FORQTS (W1,RQ1) reserves space at RQ1 for references to disk 
file Wl. 

3 CREATE 

Call CREATE (Wl, ,,,,,, ,NSECT) creates the disk file Wl with length 
NSECT sectors where a sector is 512 words. 

4. IDONE 

The statement I=ID0NE(W2) sets 1=0 if there is outstanding I/O 
to be completed to or from fale Wl, otherwise 1=1. 

5. IRANR 

The statement I=IRANR(RQ1,A(N) , NWORDS, NSECT, Wl) causes NWORDS to 
be transferred from disk file W2 starting with sector NSECT (Sector 
number 0 is the start of the file; there are 512 words per sector), 
to LCM starting at address A(N) . 

6. IRANW 

I=IRANW is analogous write statement to IRANR. 

7. FFT2 

Call FFT2(A(I1) ,B(Z2) ,N, INC) performs a fast-Fourier transform of 
length N on the real data starting at location A(I1) and the 

imaginary data starting at location B(I2) The data are mcremen- 

* 

ted by INC, -INC implies a forward transform and +INC implies a 
backward transform. 

The Mam Program 

The mam program assumes that disk files containing the first two 
Ti“l n -l n ^ Ti 

time steps, i.e., u , (u ) fc , u , and (u ) fc , exist 
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oooooo r> o o o on 


74 / 11/14 


RUN-1.CM89 


18.32.37 V5C L ARK1HA pAsfe NO. 1 


PROGRAM GAP(0UTtFSET5»FSET6»FSET7,FSET8»FSET9> 

2 IMPLICIT INTEGER (Z) 

2 LCM/BBl/DUMl (98304) 

2 LCM/BB2/DUK2 (983041 

2 LCM/BB3/DUM3 (655361 

C MAX LCM BL0CK»126976 

2 DIMENSION VlG(A 096*3*8) ,DLG (6**64**) »PLG(4096t2»4) *PZl (2048*2»24) * 

2PZ2(2''48»2i 2*) ,PZ3(2n48t2» 161 *PHAT ( Aj96,<i. 6) ,pC(64,64.6) ,WLG<*096, 

33,8) • DULG (4096,3,2) , pH AT A (*096,2,2) , DVLG (4^96,3, 8 ) 

2 DIMENSION UTNMl (24S/6) ,UTN{24576) .(JTNPI (24576) .UN(24576) , 

2VTNM1 ( 8192) f VTN ( 8192) «VTnP1 (8192) , VN(B192) 

2 DIMENSION WTNM1( 24576) tWTNU 4576), WTNP 1(24576) *WN (24576) 

2 DIMENSION USM{64,20,5) ,VSM (84,20,5) * WSM (64,20,5) , P(64*6 

?4,2) ,pZ(256,64,2 ) , Pm ( 4.96 t 2 )*PR(64,64), 0U(84,16), OV(64tl6), 

3DM ( 64 , IP) ,u(64,20,5) ,V(64,2d,5) ,W{64,20,5) ,pD(64,20,5) 

2 DIMENSION RQJ (20) ,RQ2 (20 ) ,RQ3 ( 20 > ,RQ4(20) *RQ5(20) ,RQ6(2ti) ,RQ7(20> 

2 ,RQ8 (20 ) , SPR 1 64 ) ,SPI (64) 

2 DIMENSION S(255).C(255)*S64(Z55),C64(255I 

2 DIMENSION ITAPEI2) *I05TAPE(5) ♦I06TAPEI5) *IITAP£(5) *I2TAPE(5) * 

?I3TAPE(5) ,I4TAPE(5) .I5TAPE(5) il6TAPE(5) t I7TAPE (5) * I8TAPE (5) » 

3I07TAPE(5) , I 08TAp£ ( 5 ! 

wi^ontain! 64 planes of v,v,w, each plane consists op 24 sections=i2288 

WORDS. IE' 1536 SECTIONS. PLANE 63 STARTS aT NSECT=0, PLANE I STARTS Al 
24» (1*1 ) , JF ( (24*1*1 ) .GT.1812) NSECT=24$(I*1 ) -1536 
*2 IS IOEnTICAl TO Wl , 

Dl CONSSlTS OF 6* PLANES OF DlVt EACH PLANE CONTAINS 8 SECTl0NS=4096 WORDS 
IE 512 SECTIONS. PLANE 1 STARTS AT NSECT^O, PLANE I STARTS AT 
NSECT»8*<I-1) 

PUPR CONTaINE 64 PLANES OF UPPER HaLF OF J, EACH PLaNES CONSISTS OF 8 
SECTIONS*4096 WORDS, 'PLaNE 63 STARTS aT NSeCT^o* PLANE I STARTs AT N5ECT 
=8« 1I*1)» IF( 1*1 ) *6,GT,50 4 ) NSECT=NSECT-512 
PLWR IS SAME AS PUPR 
2 COMMON OUMSM (32768) 

c EQUIVALENCE! VLGO ) »0UMl ( 1) ) * (DVLG(I) ,0UM2 ( 1 ) ) , (PL6( 1 ) ,DUM3<1) ) , 

PTdLGU ) t DUft3 (32769) ) * (WLG(l) *DOM2 (1) ) 

2 EQUIVALENCE (PZI (1) iDUMi (1) ) , (PZ2 (1 ) , OUM2 ( 1 ) ) , (PZ3(1) ,0UM3(1) ) 

2 EQUIVALENCE (DULG 1 1 > tDUMT (73729) ) , (PHAT(l) *DUM2 (I) ) , { PHAT A ( 1 ) *0UM2 

2(32769) ). (PC (1 ) »DUM2(49l5^) ) 

2 EQU I VALENCE (UTNMl ( 1 ) »OUMl (1 ) ) • (UTN H > »DUMl (24577) ) , (UTNPl (1 ) » 

2DUM1 (491531 ) (UN<1 ) »DUMl (,37 2 9) ) . ( VN(1 ) »DUMSM(1 ) ) » ( VTNMl (1 1 »0UMSM< 

34097)), (VTN(l) ♦0UMSM(el93) > * <VTNP1 (1) «DUMSM(I2289) ) 

2 EQUIVALENCE (WTNMl (1 > »DUM2 (f ) ) • (WTN t 1 » *OUM2 (24577) ) , (WTNPl Cl I * 

2DUM 2 (49153) ) , (WN(1 > tDUM2 ( 73729 ) ) 

2 EQUIVALENCE " (DUMSM(I) » SP I { 1 > ) * <Pd 

2' )>DUMSM(I )),<PZ( A «ltl ) tDUMSM(l) ) * (PA(ld > ,DUMSM (4097) 1 • ( 

3 PR (1,1), OUMSM (1 ) ) , (USM(l) ,DUMSM(4097 ) ) , (VSMfl ) .DUMSM (10*97) ) * 

4 ( WSM ( 1 ) » DUMSM ( 1 6897) ) » ( OD ( 1 , 1 ). DUMSM ( 1)>»( DV ( 1 , 1 ), DUMSM 

5(1(123)), ( 0W{1,1) *DUMSM( 2049) ) , <(J ( 1 t 1 > 1 1 »DUMSM ( 3073) ) , ( V ( 1 • 1 » 1 > 

6»DUMSM ( 9473 ) ) , (rf( 1,1,1) , DUMSM ( 1 5873) ) . (PD ( 1 * 1 * X ) .DUMSM (22273) ) 

DATA RQ l/20*g./» RQ2/2Q*Q./* RQ3/2n*n ./ t RQ4/2o*0 •/• R05/20" 0 ./ 

DATA RQ6/'20*0,/’, RoT/EOaO./.RQS/ZOttS./ 

DATA Wl ,W2,Dl .PUPR.PLWR/fLWl ,2LW2,2LDl ,4 lPUPR,4 l p u WR/ 

DATA w3,W4,W8/2LW3»2LW4,2lWB/ 

DATA ITAPE/1.4LT4PE/* I05TAPE/**3«0,8lXX014477/,i06TAPE/*»3 #0 * 
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RUN-LCM89 


GAO 


74/11/14 


16.32.37 


VBCl_ARKlHA PAGE NO. 2 


2 

4 

6 

10 

12 

14 

16 

2fl 

22 

24 

31 

34 

46 

61 

73 

106 

120 

133 

145 

160 

172 

205 

217 

232 

2A4 

257 

271 

304 

332 

334 

337 

341 

343 

347 

351 

354 

361 

362 
371 
371 
400 
402 
404 

406 

407 
422 


44 0 


28LXX01 A883/, IlTAPE/A*3*0,8LXX005ll2/, l2TAPe/4*3*0|8uXX00576q/ 
DATA 13TAPE/4,3#J,8LXX005663/.H7APE/4.3*0.8LXX00532?/,I5lAPE/4, 

2JO0*BUXX00703<)/* I6TAPE/4*3PO«8LXX009091/*17TAPE/5 # 0/ 

3* I8TAPE/5»0/ 

DATA I07TAPE/5*0/, IO8TAPE/5*0/ 

CALL MEMREQ(262144.11 
CALL F0RQT5(Wl»R3ll 
Call forotS(*2*r32> 

CALL F0RQTS(W3,R33) 

CALL F0RQTS(K4,R34) 

CALL F0RQTS(W8,R38) 

CALL FORQTSIqI »R35) 

CALL FOROTS ( PUPR*RQ6) 

CALL F0RQTS[PlWR,RQ7> 

00 3 Ini, 8 
3 DUM1 ( I ) =0 . 

CALL CREATE (HI ,U*HT*0*0*0»0*0*0,1536> 

I»lRANI«((RQl»oUHl (1) • 1* 153f »W1) 

CALL CREATE (W2 ,UfRT*0|0,0»0«0, 0,1536) 

I=IRANW(Rq2,oUH1(2) ,1*153=, W2) 

Call create(«3 ,u*rt.o,o,o,o,o,o,is36) 

T=IRANW(Rq3,OUh1 (3) ,1 ,153^,W3) 

CALL CREATE (W4 ,U* RT, 0,0, 0*0*0* 0,1536) 

I=JRANW(RQ^,DUMl (A) ,l,153=,w4) 

Call creatE<w8 ,u*rt*o*o,o,o,o,o,i536) 

I=IRANW(RqS,DUh1 (8) , 1*1535, W8) 

CALL CREAt£(D1 ,U*HT*0, 0,0, 0,0, 0,512) 

l=IRANW(RQS,OUMl (5) ,1,511 , dm 

Call create (pupr,u*rt*o,o,o, 1 , 0 , 0 , 512 ) 

TsiranW(RQ6.0Uh 1 (6) , 1*511, pUPR) 

CALL CREATE (PLWR,U*RT*0,0*0,0,i,0,sl2j 
I*1RANW (Rq7 ,011m 1 (7) 1 1 *511 ,PLHR) 

912 1*IdONE (WU *IDONE ) *IDONe (W3 ) * IqONE (WA ) *IdONE tW8 I ♦ IDONE (D 1 > 

= MDON£<PUpR)*lDONE(PLWR) 

IF ( l»NE.8) GO TO 912 

CALL 0PEN(5lFSET5i0. 2340008) 

00 920 kb 1,8 
NSECT c (K-l)4lg2 - 
M=MOLHK*2) 

00 9li J=1 *3 

JK=(U-1) <*32768 + 1 

ftEAD<5) (DUMSM(I) *I»1 *32768) 

1F<M.E0,0)G0 TO 9°» 

Small* out(0umsm(1) Jdumi (jk> *32768) 

50 TO 910 

9 j9 Small out<oumsm<i> *dum2(jk> * 32768) 
glo CONTINUE 
911 I=ID0NE(W1) 

IF(I.NE.l)SO TO ?11 
IF (M.EQ.O ) GO TO 919 
I=IRANW(RQl*DUMl (1) *98304, NSECTtWl ) 

GO TO 9^0 

919 I=IRAN*(RQl,DUM2(l> *98304, NSECTtWl) 

920 CONTINUE 

Ifirst»o 


„, r n v U 3 & 


^PEODUCIBILITY C" 

Pit " *3 W 
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441 

443 

444 
446 
452 
454 
457 

464 

465 
474 
474 
503 
505 
512 

514 

515 
530 
530 

544 

546 

5S0 

553 

555 

557 

563 

565 

570 

575 

576 
605 
605 
614 
616 
623 

625 

626 
641 
641 
655 
657 
661 
663 
665 
671 
673 
676 

703 

704 
713 
713 
722 
724 
731 

733 

734 
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IF(IFIR$T. EQ. 1 ) GO TO 941 
00 940 K«l ,8 
NSECT=(K-1')*192 
M=*M0D(K.2) 

DO 93, J=l,3 

JKa< J-1)*3276B*1 

REA0IS>1 (DUMSMU) .1*1*32768) 

1FIH.EQ.OJGO TO 929 

SMALL OUT (GUMSM ( 1 ) *DUM1 (JK) .32768) 

GO TO 930 

929 SMALL OUT(OUMSMti)»DUM2UK)»32768) 

930 CONTINUE 

931 I=«ID0NElW2>*ID0NE(Wl> 

IF ( 1.NE.2) GO TO 931 
IF(M.£Q.O)Go TO 939 
jaIRANW(RQ2tDUHl (1) . 98304 , NSECT. W2 ) 

SO TO 940 

939 I® lRANhMRQ2» DUM2 ( 1 ) . 98304. NSECT, W 2 > 

940 Continue 

941 CALL AFSREL(5LFSET5) 

CALL OPENI5LFSET6.0. 2340008) 

OO 960 K=l »8 

NSECT=(K-1)*192 

M=MOl)lK.2) 

DO 9§. J=l*3 

JKo(J-l)*32768*3 

fiEADlb) (DUMSM(I) .I«t .32768) 

lF(M.EQ.0)GO TO 949 

SMALL OUT (DuMSM 1 1 ) »QUM1 ( JK) .32768) 

GO TO 950 

949 SMALL OUT(OUMSM(l ) .0UM2tJK) *32768) 

95 0 CONTINUE 

95] I=ID0NE(W3)*ID0NE(W2) 

IF ( I .NE.2 ) GO TO 951 
IF (M.EQ.O >GO TO 959 
I3IRANW(RQ3.DUM1 (1 ) ,98304, NSECT iW3) 

GO TO 9?0 

959 I^IRANW (RQ3.DUM2 ( 1 ) .98304. NSECT.W3) 

960 CONTINUE 
IFIIFIRST.EQ. 1 ) GO' TO 981 
DO 980 K=1 .8 
NSECT»(K-1)«192 
M»M0U(K, 2 ) 

00 9?; J=1 .3 

JKa (J-l)*32768*l 

READ (6) (DUMSM(I) .1*) .32768) 

IF I M.EQ.O ) GO TO 964 

SMALL OUT (DUMSM (1 ) .IJUMI (JK) .32768) 

GO TO 970 

969 SMALL OUT(DUMSMU) .0UM2(JK). 32768) 

97n CONTINUE 

97] I=I00NE(W4)+I0QNE(W3) 

IF ( I .NE,2) GO TO 971 
IF (M.EQ.O)GO TO 979 
I3IRANW1RQ4,0UM1 (1) .98304, NSECT.W4) 
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747 


GO TO 980 

747 

979 

Is I RANH {RQ4,DUM2 ( 1) ,98304 t NSECT,M4) 

763 

9Bo 

CONTINUE 

765 

9B 1 

CALL AESREU5LESET6) 

767 


DO 1 1=1,32768 

774 

1 

oumsm(D=o. 

776 


W=3, l4lb926535698/2o4a. 

777 


UO 4 i=i,255 

1001 



1 0 1 0 

4 

C(D*COS ( (I»l>**f) 

1021 


fl=3. 1415926555898/32. 

1022 


UO 5 1=1,255 

1024 


564 ( i ) =SIN ( ( I“1 1 »W) 

1033 

5 

C64(I)=C0S( !I-1)»WI 

1045 


0ELT=.0073 

1046 


0TDT»UELT/l2. 

1047 


TTDT»2.*DELT/3. 

1051 


DELZ=20,/64. 

1053 


C 1 = 1 . / ( 12.*DELZ) 

1055 


X=2r, 

1057 


66=7456540. 4*2/X**2 

1061 


C r 2aCl 

1062 


CN*. 1453 

1064 


C7=12./<5.*DELT) 

1066 


C8«l,/(12.*OELZ> 

1071 


Cll=l./(12.#0ELZ»*2t 

1073 


TDT=2,poELT 

1074 


ITIME=4o 

1075 


iset=_ 

1077 

1001 

CONTINUE 

1077 


Ih=H0D(ITIME.4141 

1104 


IMITIHE.LT, J)GO TO 400 

1105 


IF(MOO<ITlMe,e) , 3T.3) JSET»1 

1113 


IF ( I SET.EQ, 0) GO To 4i;0 

1114 


00 873 i=l»6i 

1121 


SPR (11=0. 

1122 

873 

SPI ( I) =0. 

1123 


RMSD=.. 

1123 


RMSU*0 . 

1124 


RMSV»0« 

1125 


RMSW=o» 

1125 


SK1=; . 

1126 


SK2=0 . 

1126 


umax=u. 

1127 


VMaX»0, 

1127 


WMAX»0. 

1130 


CRMSUsO. 

1130 


CRMSVsO, 

1131 


CRMSrf=o . 

1131 


CUMAXoO, 

1132 


CVHAX=0. 

1132 


CWHAX=0 , 


C 

C 

»«„*«..**»* PHASE I <h»***»4o«n***44 


C 
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1133 

1161 
1163 
1166 
1166 
1167 
1172 
1205 
1207 
1217 
1232 
1246 
1246 
1261 
1275 
1275 
1310 
1324 
1324 
1337 
1353 
1366 
137 0 
1372 
1401 
1410 
1415 
1420 
1430 
1440 
1454 
1454 
1470 
1470 
1504 
1504 
1520 
1522 
1535 
1537 
1541 
1543 
1547 
1547 
1551 
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C START WITH PLANES 63,64,1,2,3*4,5,6 IN LCM 

c v 1 1 ) at segment », v<2) at segment 24* vin> at segment 24 *n«ij 

c VC64I at segment 1512, EiC H V CONTAINS a 22q8 WORDS 


c 

VARIABLE LOCATIONS as 

follows at 

start of 

c 

TIME* 1 

2 

3 

4 

c 

W1 JtN-1) 

UTtN) 

UT(N-l) 

U(N) 

c 

W2 UTtN-1) 

UIN) 

utN^n 

UT<N) 

c 

W3 U (N) 

U(N*1> 

UT(N) 

UT(N-l) 

c 

W4 UT (N) 

UT ( N- 1 > 

UIN) 

UtN-1) 


400 IcIDDNE(Wl) ♦ID0NE1W2) *IDONE(W3) *IDONE(W4)+IDONE<W8)+IOONE{D1) 
?*IDONE (PUPR) *1D0NE<PLWR) 

IE { I ,NE,8) GO TO 400 
IF(ITIMt,LT,0)GO TO 1901 
IE(ISeT;EQ.O!GO TO 1902 
00 155 N C 1 *8 
NSECT*(N-1)R192 

I3i> I»I00NE(W1)+ID0NE(W2)*ID0NE(W3 ) *IdONE<W4) 

IF(I.NE,4)G0 TO 130 
GO TO <131,132,133*134) IM 
131 I“IRANR(RQl f VL0<l)*983J?4*NSECT«Wl) 

IbIRANR<RQ4»dVLG(1) *903 0 a *nSeCT *W4) 

GO TO 139 

! 32 IsIRANR(RQ3,VLG(i) *9B304*NSECT,W 3 ) 

I=IRANR(RQ1* dVLGI D . 963oJ.NSECT.Wl) 

GO TO 139 

133 I=IRANR<RQ2,VLG<1) *98304. N$£CT*W2) 

IbIRANR(RQ3,0VLG<1) *9830;*NSECT*W3) 

GO TO 139 

)34 1=IRANR{RQ4* VLG ( 1 ) * R8304 ,NSECT*W4 ) 

I*IRANR (RQ2, DVLG < 1 > ♦ 9830T * NSECT * W2) 

139 I=IOONE(W1)*IdONE{W‘ : )*IDONe<W3)*IoONE<W4) 

IF ( I ,NE,4) GO TO 139 

00 141 K=l,a 

SMALL IN ( DUMSM ( 1 ) ,VLG(1.1,K) *12288) 

SMALL IN<DUMSM(13289) ,OVLS< *1*K)*12288) 

00 140 1=1 *1^288 

140 DUMSM < 1 1 =DUMSM < I ) +T0T*0UMSM (1*12288) 

141 SMALL OUT (DUMSM(l) .VLGd.l.K) *12288) 

GO T0(151*IS2«153.ib4) 1M 

151 1 = IRANW<RQ1,VLG(1> *983o4,NSECT*W1 ) 

GO TO 135 

152 I=IRANH(RQ3,VLG<1) *98304. NS£CT»W3) 

GO TO 1R5 

153 I«IRANW{RQ2,VLG(1) *98304. NSECT.W2) 

GO TO I 35 

]54 I s IR ANW 1 RQ* * VLG ( 1 ) . 98304 *N5ECT*W4) 

155 CONTINUE 

156 r»IDONE(Wi)*IDONE<W2)*IOONE (W 3 ) *IOON£(W4) 

IF ( I ,NE,4) GO TO 156 

157 I^IDONEtOlT 

IF ( I »NE, 1 ) GO TO 157 

IF(ITIME,GT.3)G0 TO 600 

ID«1 

IMPoIM 

GO TO 499 
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1552 1911 IMP=4 

1553 10=1 

1554 GO TO 499 

1555 1902 XFariHE.GT.3)G0 TO 600 

1561 IF(IH.EQ.l) IMP'4 

1564 IFUM.EQ.2) IMP»l 

1567 IF<IH,EQ.3J IMP=2 

1572 IF<IM.£G.4> IKP=3 

1575 10=1 

1576 499 CONTINUE 

1576 IF ( ILJ.NE.l ) I=IRANR (RQ5»DUM3 < 1 ) »32760»O»Ol) 

1613 00 T6(*61,*6f»* 6 3,464) IMP 

1623 461 I = IRAnR(Rq 3,0UM1 U> *9S^CI4»rnW3) 

1637 ' GO TO 465 

1637 462 I=IRANR(R02,OUHl (i) ,98304,0*^21 

1653 GO TO 4»5 

1653 463 I=IRANR(RQ4,DUM1(1) ,98304, n,W4> 

1667 GO TO 465 

1667 464 UIRANR(RQ1,DUM1 (l) ,98304, o»Wll 

1703 465 I=IDONE(Dl>,IOONE(WA|+IOONE(W2)*rDONE<W3)*lDONE(W4> 

1721 IF<I.NE.5)G0 TO 465 

1723 00 490 M=l» 1 6 

1724 I>{m6d(H,2) .EO.OIGO TO 480 

1730 NSECT=M«96 

1731 NWOROS*98304 

1733 J=1 

1734 IF(M,EQ,151NW0R0S=49152 

1737 475 I»ID0NE(W1)+ID0NE(W2) ♦ I DONE (W3)*ID0NE(R4> 

1752 IF(I.NE.4)G0 TO 475 

1754 t 00 TO (291,492*493*4941 IMP 

1764 *491 I*IRANR(RQ3*DUM2(J) ,NW0RDS,NSECT,W3) 

2000 80 TO *^ 3 

2000 492 I=IRANR(RQ2,0UM2( J) ,NWORDS,NSECT,W2) 

2014 50 TO *? 3 

2014 493 I=IRANR(RQ4,DUH2(JJ ,NH0R0S,NSECT,W4) 

2030 50 TO *73 

2030 494 I=IRANR(RQl,DUM2(J),NW0RDS,NSECT,Hl) 

2044 473 IFMM.NE. 15), OR. INSECT. EQ»"))GO TO *79 

2053 NSECTsO 

2053 J=49i53 

2054 GO TO 475 

2055 479 I=ID0N£(Dl) 

2057 IF ( I ,NE. 1 ) GO TO 479 

2061 IFIM.EQ.llGO TO 477 

2063 J=M*i 

2063 NSECT=(J/2-2)P64 

2067 IOA=l 

2070 IF(M0D(J,4l«NE.ij) 104=32769 

2075 I*IRANW(RQ5,0UM3( IOA), 32768, NSECT.ol) 

2111 GO TO 4(7 

2111 48o NSECT»M*96 

2113 NWORDSP98304 

2115 476 IpIO0NE(W1)*IOONE(W2)*IOONE(W3)+IOONE(W4) 

2130 IF(I,n£.4)G0 TO 476 

2132 1F(M,£Q« 1ft) GO TO 4T7 


<\/ 

xJ> 
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2134 

2144 

2160 

2160 

2174 

2174 

2210 

2210 

2224 

2226 

2230 

2232 

2234 

2234 

2242 

2245 

2261 

2261 

2277 

2301 

2303 

2305 

2320 

2336 

2340 

2341 

2342 

2351 

2352 
2355 
2360 
2363 

2366 

2367 
2370 
2373 
2376 
24 01 

2404 

2405 
2405 
2405 
2407 
2411 

2437 

2441 

2442 
24 52 
2462 
2476 
2476 
2512 
2512 
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SO T0<495, 496, 497,496) IMP 

495 UIRANR(RQ3,DUM1 (1) ,NWORDS, NSECT, W3) 

SO TO 478 

496 l“IRANfi(RQ2,0UMl(l) ,NM0R0S.MSECT.W2) 

GO TO 478 

497 I»IRANR_(RQ4,0UM1<1) ,NWORDS, NSECT «W4> 

GO TO 478 

498 IcIRANRIRQl ,DUM1<1) ,NWORDS, NSECT, Wl ) 

478 I=ID0NE(D1) 

IFtl.NE.DQO TO 478 
IF ( IU, EG, 1 rGO TO 477 
IF<M;fcG.l6)G0 TO 477 
IOA°l 

IF <MOD (M»4) »NE, 0) I0A»32769 
NSECT=(M/2)*64 

I=IRANR(Rq 5,DUM3(I0A) , 32766 , NSECT, Dl ) 

477 CONTINUE 

CALL DVRGNC tCl2,RMSD»RMSU,RMSV,RMSW.SKl.SK2,UMAXtVMAX,WMAX,M,IO,OE 
?LT> 

49(i CONTINUE 
484 I*IDONE(Ol) 

IF < I.NE.l ) Go TO 484 

I=IRANW|RQ5,DUM3 (32769) ,32768,448,01) 

489 I = IoOI9E(W1)*IOONe(W2)*IDONe1W3)*IDONE<W4)*IoONIE(D1) 

IF ( I.NE.5) GO TO 489 
IF ( ITIHE.LT ,0) SO TO 601 
IF I ISETiEQ. o)SO TO 608 
GO T0(501,502,503) ID 
501 IDs2 

IF(IH.EQ.l) IHP=2 
IF (IM.EQ.2) IMP=3 
IF ( IM, EQ.3) IMP»4 
IF (IM.EQ.4I IHP*1 
GO TO 499 
5fi2 IO»3 

IFtIH.EQ.l ) IMP*3 
IF ( IM, Eq.H) IHP n4 
IF (IH.EQ.3) IMPcl 
IF(IM, EG, 4) IMPoa 
60 TO 499 
503 CONTINUE 

loo Continue 

46n lalOONE(Ol) 

IF (I.NE.l)GO TO 460 

600 IsIDONE (Wl ) +IOONE (W£) ♦ I DONE (W3) +I 0 ONE (W4) *I00NE (W8) *100N£(Dl) 

2* I DONE (pUpR) *IDONE(PLWR) 

IF ( I «NE. 8 ) GO TO 60J 
IF < iSET.EQ. 0 ) GO TO 0 O 8 
GO TO ( 6 1, §02, 603, 604) IM 
608 GO TO (602,603,604,601) IM 

601 l*IRANR(RQl,VLG(l) » 983^4 , o,W1) 

GO TO 605 

6(i2 I = IRANR(RQ3»VLGU) *983o4,o.W3) 

GO TO 605 

603 I*IRANR(RQ2,VLG(1) ,98304, o,W2) 
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2526 GO TO 605 

2526 604 I=1PANR <RQ4, VLG C 1 ) ,98304,o,W4) 

2542 605 I=IOONE(W1)*IDONE(W2)+IDONE(W3>*IDONE(W4) 

2555 IP(I.NE,4)60 TO 605 

2557 ZH2=1 

2560 lOAcl 

2561 00 620 IPLANEb1.64 

C tRANSFER ROWS 45 TO 64 TO CALCULATE 47 TO 62 

2562 N2=28l7 

2563 N3=1280 

2563 N4=1 

2564 NEXT°1 

2565 M»1 

2566 613 K«=ZM2-1 

2570 00 621 Z*l»5 

2572 K»K*1 

2573 IF(K.0T.S)Kal 

2576 SMALL INI U ( 1 1 N4*Z) »VLG(N2»1»K) tN3) 

2612 SMALL INI Vtl*N4,Z! «VLG(N2.2.K) «N3) 

2626 621 SMALL INt W(1 ,N4,Z) ,VLG (N?«3»K) ,N3) 

2644 60 T0<622.62J,624»623)NEXT 

2654 622 CALL UVOTmP (U. V • W.OU.DV , DW ,M,CN, C 0, Cll • Cl2* lTIME»RMSO*RMSU»RMSV , 

2 RMS W* SKI *SK2*UMAX*VmaX»WhAX*SPR* ISET) 

2704 SMALL OUT <DU t 1) . DVLG(29*S;i, JOA) t 1 Q24) 

2714 SMALL OUTlOVIll ,DVLG(2?45,?, IOA) ,1024) 

2723 SMALL OUT(Dwll) *DVLG(2945»3,I0A) ,1024} 

2732 DO 614 j=l,4 

2733 J J=J* 16 

2734 Do 614 1=1,64 

2736 DO 614 Z«=l,5 

2752 U(I»J»Z)=U(I»JJ»Z> 

2754 V(I,J,Z)=VU,JJ,ZJ 

2755 614 Wl'lt j»Zl=W<I»JJtZ) 

c calculate rows 63,64 and i to i* 

2763 N2=l 

2763 N3=l 024 

2764 N4=5 

2765 NEXT=2 

2766 M a 2 

2767 60 TO 613 

27To 623 CALL DVDTMP |U* V, rf, DU,DV.DW,M,CN»C8»Ci 1 *C12* ITIM£,RMSD»RMSU»RMSV, 

?RMSVJ*SKl*5KH,UMAX,VMAX,tfMAX»SPR»lSETl 
3020 SMALL OUTlDU(l) ,DVLQ(3969 t l ,I0A1 ,128) 

3030 SMALL OUT (DV tl) «DVLet3969»2,IOA) ,128) 

3037 SMALL OUT <DW ( 1 ), OVLG (3969,3 , 10A) , 1 28) 

3046 SMALL OUT (OU < 129 ) , DVL6 ( 1 , 1 , 10A ) ,896 1 

3055 SMALL OUT ( DV ( 129) ,DVLG ( 1 ,2 , IOA) ,8961 

3064 SMALL OUT ( DW 1129) * DVLG I 1 ,3 , IOA) ,896) 

C CALCULATE ROWS 15 TO 3JJ 

3073 N2=769 

3073 N3=»lHeO 

3074 N4=l 

3075 NEXT»3 

3076 * M=3 

3100 60 TO 613 
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3100 

3130 

3140 

3147 

3156 

3156 

3157 

3160 

3161 
3163 
3163 

3213 

3223 

3232 

3241 

3244 

3246 

3264 

3266 

3270 

3271 

3272 
3275 
3300 
3304 

3306 

3307 
3317 
3327 
3344 
3344 
3361 
3361 
3376 
3376 
3413 
3415 
3417 
3423 
3440 
3442 

3445 

3446 
3451 
3453 
3471 
3473 
3511 
3525 
3527 

3531 

3532 
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624 CALL DVOTMP (U, V , W,OU,DV, DW, M,CN, C8, Cl 1 ,C12. ITIME , RMSD.RMSU.RMSV, 
2RMSW,SK1 »SKZ,UMAX.VMAX»WmAX.SPR» I SET) 

SMALL OOHDU(l) .DVL6(897,I.I0A| .1024) 

SHALL OUT (DV<1)»DVL0(89?,2,ioA) ,1024) 

SMALL OUT (DW (1 ) , DVLGI897, 3. IOA) .1024) 

C CALCULATE ROWS 31 To 46 
N2*l 793 
N3“1280 
N4=l 
NEXT*4 
M»4 

00 TO 613 

625 CALL DVOTMP (U. V . W. DU.OV. DW. M, CN.C8.CU . C12. I TIME .RMSO.RHSU.RMStf , 

2RHSW. SKI. SK2. UMAX .VmaX»WmAX.SPR» ISeT) 

SHALL OUT ( DU ( 1) , DVLG ( 1921 ♦ 1 . IOA1 , 1 02*> 

SHALL OUT1DVU) » DVLG I 1921 .2 » JOA ) ,1024) 

SHALL OUT(DW(l) ,DVL8(1921.3.I0A) .1024) 

Io MODtIPLANt.,2) 

IF(I.NE.0)GO“TO 619 

6C6 I=I00NE(W1) +IDONEIW2) *1D0NE<W3) *ID0NE<W4) *IOONE (W8) 

XF(I.NE,5)60 TO 606 

IF(IPLANE.EQ,64) GO TO 616 

IO0=ZH2-1 

I=IPLANE*5 

IF (I.6T.64) IeJ-64 

NSECT»2?* ( I ♦ 1 ) 

IF INSECT. GT. 1512) NSECTaNSEC T -1536 
IFIITIME.LT.OIGO TO 631 
IFIISET.EQ.OIGO TO 630 
GO T0I631 ,632,633,634) IH 
63r> GO TO (632,633,634,631 HH 
631 I=IRANR {RQlivLGtl.l.IOB) , 24576, NS ECT ,W1 ) 

GO TO 616 

63g I=IRANR {R03,VLG(l»i»I08) .24576.NSECT.W3) 

GO TO 616 

633 I=IRANR <RQ2,VLG(l.),IOB) .24576.NSECT.W2) ' 

GO TO 616 

634 I=IRANR (ROa.VLG(i»i.IOB) .24576.NSECT.WA) 

6 1 6 I08=I0A-1 

NSECT=IPLANE*24 

IFtNSECT.GT,1512)NSECT=NSECT-1536 
I=IRANW ( RQ8 , DVLG ( 1 . I • IOB) , 24576 .NSECT.W8) 

61 9 IOAbJOA.1 “ 

IF (lOA.GE.S) IOA=1 

ZH2aZM2*l 

IF<ZM2,GE.9)ZM2»1 

620 CONTINUE 

626 r»IDONE(Wi) *IDONE(W2) *IDONE <W3) +ID0NEtW4) ♦ I DONE (W8) 

IF (I *NE,5) GO TO 626 

rFUTlHE.LE.3)I>=IRANRCRQ5,Di-G<U ,16384.0,01) 

. I = IRANR'(RQ8,vlG(1) ,98304, 0.W8) 

1'9 I=ID0NE(W8) 

IF ( I «NE. 1 ) GO TO 1J9 
ZH2=l * 

IOA®! 
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3533 

353* 

3537 

3546 

3546 

3553 

3555 

3555 

3556 

3557 

3560 

3561 

3562 
3564\ 

3566 

3567 
3572 
3606 
3622 
3640 
3656 
3653 
3663 
3666 

3670 

3671 
3673 
3707 

3711 

3712 

3720 

3720 

3721 

3722 

3723 

3724 

3725 
3730 
3740 

3743 

3744 

3745 

3746 

3747 
3750 
3750 
3753 
3763 


3766 

3767 
3770 
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00 20 IPLANE*1 *6 * 

C TRANSFER ROWS 45 TO 64 TO CALCULATE 47 TO 62 

IF { I TIME. GT .3) GO TO 26 
SMALL INtP(l) .DLGll, l»IOA) .4096) 

30 TO 28 

26 00 Z< 1=1.4096 

27 P < I J =0. 

28 CONTINUE 
N2=28l7 
n 3«1260 
N4=l 
NEXT*1 
Mel 

13 K=ZM2-1 

00 21 Z=1,S 
KeK+1 

IF(K.GT.8)K=1 

small IN ( USM ( 1 « N4 « Z ) * VLG [NHt 1*KJ ,N3) 

SMALL IN(VSM(1,N4»Z) .VLG(N2»2*K> .N3) 

2t SMALL IN(WSM(1,N4,Z) .VLG(N2*3*K) ,N31 
GO T0I22.23* c 4*25)NbXT 

22 lF(ISET.EQ.O)C7*l./TOT 

CALL CALCPR (USMi VSMi WSMip(M»Cl 1 C7 > 

C7=.12./(5.AUELT) 

00 14 jal ,4"' 

JJ=J*16 
00 1 * 1=1.64 
DO 14 Z=1.5 
USMU»J»Z)=USM(I,JJ.Z) 

VSM(I.J*Z)=VSM(I»JJ.Z) 

14 WSMd.J.ZleWSMU.JJ.Z) 

C CALCULATE ROWS 63.6a AND 1 TO 1* 
nHbI 
N3=l 024 
N4o5 
NEXT*2 
M=2 

GO TO 13 

23 lF(IS£T,EQ.o>C7»l./TDT 

CALL CALCPR(USM.VSM.WSM,p,M»Cl.C7) 

27=12,/ (5.*UELT) 

C CALCULATE RoSs 15 TO 30 
N2=769 
N3=l280 
N4= l 
NEXT =3 
M=3 

GO TO 13 

24 IF (ISET ,EQ,0)C7=1 ,/TDT 

CALL CALCPR(USM»tfSM»WSM»p,M»Cl*C7) 

C7=12./(5.*UELT) 
c calculate rows 31 to 46 

N2=l793 

N3=1280 

N4=l 
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3771 NEXT-4 

3772 M=4 

3773 30 TO 13 

3773 25 lF( ISET.EQ, o>C7=1,/TDT 

3776 CALL CALCPR (USM» VSM,WSM,p,M* C1»C71 

4006 Q7-12,/ (5, to 0£LT) 

4011 CALL FFTFXYlP) 

4013 I-MOO (I0A»2) 

4017 IF(IiN£.0)G0 TO 16 

4020 I0b-IOA-1 

4022 SMALL 0UT(Pd*l *1),PLG(1 *2* IOB) *2048) 

4031 SMALL 0UT(P(1*1 ,2) .PLG(2G49,2, IOB) ,2048) 

4040 SMALL 0UT(P(1,33,1) ,PLS(1 ,2, IOA) *2048) 

4047 SMALL OUT (P (1 *33*2) .PLG (2049,2, J04> .2048) 

4056 GO TO 17 

4u56 16 IOB-IOA+1 

4060 SMALL 0UT(P(1,1 ,1),PLS(1 , 1 , IOA) ,2048) 

4067 SMALL OUT(P(l,l ,2) , PLG (2049,1 ,IOA) ,2048) 

4076 SMALL OUTtPll, 33,1), PLGIl *1 . IOB) *2048) 

4105 SMALL 0UT(P(1,33,2) *PLG(2049,1,IOB) *2048) 

4114 17 CONTINUE 

C PAUSE IF THERE IS ANY OUTSTANDING I/O 

4114 I=MOD(IPLANE,2) 

4120 IF ( I »NE* 0 ) GO TO 19 

4121 113 I -I DONE (01 T *IDON£.(WB) 

4126 IF { I «N£,2) GO TO 113 

4130 IP (IPLANE.EQ, 64) GO TO 116 

4132 I=IPLANE*3 

4133 IF(I.GT. 64)1=1-64 

4136 NSECT=8»(I-1) 

C ?EAO OIV OF PLANE+2 INTO LCM OlV(l,l,IOA) 

4140 IOB=IOA-l 

4142 TF (ITIME,LE,3) I- IRANR (ROg* DLG (1*1*108) *8192,NSECT , Dl ) 

4162 I=IPLANE*3 

4164 IF’<I,GT,64) Id-6 4 

4167 NSECTs=24*(I*1) 

4172 IF(NSECT.GT.1512)nSecT=NSecT-1S36 

C HEAD VELOCITIES OF PLANE*4 INTO LCM VEL (4096* 1 *ZM2> 

4176 IOB-ZM2-1 

4200 tF(I0B.LE,O) STOP 

4203 I = IRANR(RQ8,VLGd, 1*108) ,24S76*NSECT,W8) 

4220 116 I=IDONE(PUPR) 

4222 IF(I.NE,l)GO TO 116 

4224 117 I-IOONE(PLWR) 

4226 IF ( I • N£, 1 ) GO TO 117 

4230 NSECT-84(IPLANE'2) 

4232 TF INSECT. GT, 504 ) NSECT=NSECT**512 

C */RlTE(PLG(l*l,IOA) FROM LcM TO DISK 

4236 IOBoIOA-1 

4240 I=IRANW (RQ6*PLG (1*1 ,IOB) «8192*NSECT*PUPR) 

4253 I-IRANW (RQ7*PL6(1*1*I0A) ,Bl92,NSECT,PLWR) 

4270 118 CONTINUE 

4270 19 10A=I0A*1 

4272 IF (I0A,GE.5)I0A=1 

4275 ZM2bZM2*1 
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4276 

4301 


4303 

4331 

4333 

4346 

4350 

4352 

4365 

4367 

437) 

4404 
4406 
44l 0 

4411 

4412 

4414 

4417 

4421 

4422 
4433 
4443 
4454 

4466 

4467 
4470 
4501 
4512 
4522 

4524 

4525 
4536 
4546 
4557 

4571 

4572 

4573 
4604 
4615 
4621 

4623 

4636 

4640 

4642 

4655 

4657 

4661 


GAB 

IF (ZM2,GE»9)ZM2=1 
20 CONTINUE 

0 PHASE II et.H.toootoOi). 

C WAIT FOR OUTSTANDING I/O 

1 49 I = ID0NE ( PUPR) .IDONEIPI-WR) * I DONE < 0 1 > + I DONE ! Wj ) ♦IQONE (W2 1 + I DONE t W3) ♦ 
b I DONE (W4)*ID0NE(W8! 

IF ( I «NE.9) GO TO D* 

C READ INTO PZ1.PZ2.PZ3* PRESSURES FOR I»l» 64, J= 1 , 32»Z=1 . 64 
I=IRANR (RQ6.PZ! (1.1.1) .98304 .0 ,PUPR) 

26) I»IDONE (PUPR) 

IFU.NE.DGO TO 261 

I=IRANR (RQ6«PZ2<1»1»1) .98304 .192.PUPR) 

262 I»lDONE(PUPRl 
IFU.NE.D60 TO 262 

I=IRANR (RQ6.PZ3U.1.1) .65536 .384.PUPR) 

C WAIT FOR OUTSTANDING I/O 

263 I*IDONE(PUPR) 

IFU.NE.1JGO TO 263 
MM=0 

36 M=1 

UO 30 J«1.8 
JK* ( J« 1 ) *256*1 
UO 32 K-1.2A 
KK=K*24 

Small in(pzu.k .u .pzdjk.i.ki .256) 

SMALL INIPZU.K .2) .PZKJK.2.K) .256) 

SMALL IN (PZ 1 l.KK. 1) .PZ2(JK, l.K) .256) 

32 SMALL IN(PZU.KK»2>.PZ2(Jk*2*K)i256) 

DO 33 Kb 1 . 1 6 

KK=K‘48 

SMALL IN(PZU.KK.l) . PZ3 ( JK. 1 . K) .256) 

33 SMALL IN(PZ(1.KK.2) »PZ3(jK*2»k) .256) 

CALL DELSQ (PZ.C.C64. S.S64.M.MM.C6) 

DO 34 Kal.24 

KK=K*24 

SMALL OUT (PZ (1 .K .1) .PZ1 (JK.l.K) .256) 

SMALL OUT (PZ (l.K ,2) .PZl ( JK.2.K) .256) 

SMALL OuTfPzU.K^.l) .PZ2{JK,1,K) ,256) 

34 SMALL OUT (PZ (1.KK.2I .PZ2 (JK.2.K) ,256) 

DO 39 K«1.16 

KK=K*4S 

SMALL OUT(PZ(l.KK,l) .PZ3IJKU.K) .256) 

35 SMALL 0UT(PZ(1.KK.2) .PZ3(jK.2,K) .256) 

3b M»H»* 

IF(MM.EQ. 32)60 TO 38 

C WRITE PZI.PZ2.PZ3 To PRESSURES FOP 1=1 .64. J=1 .32.Z*1 . 64 
I=IRANW (RQ6,pzl (1.1.1) .98304 ,0 ,PUPR) 

161 I=IOONE (PUPR) 

IF ( I «NE, 1 ) 60 TO 161 

I=IRANW (RQ6.PZ2I 1. 1*1) .98304 .192. PUPR) 

162 I=IDONE(PUPR) 

IFU.NE.D60 TO 162 

I = IRANW (RQ6.PZ3(D1 .1) *65536 .384. PUPR) 
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C 


467A 

4676 

C 

163 

4700 

4713 

4715 

4717 


171 

4732 


172 

4734 

4736 

C 


4751 

4753 

4755 


173 

4756 

4757 

C 

38 

4757 

4773 

4775 


181 

4777 

5012 


182 

50 1 4 



5016 

c 


5031 


183 

5033 

5035 

5053 

c 

191 


c 



c 


5055 



5055 

5056 


150 

5104 

c 



c 



c 


5106 

5121 


256 

5123 

5125 



5126 

5130 


203 

5132 

5134 

5136 


204 

5140 

5154 



5170 



5172 

5175 


205 


WAIT FOR OUTSTANDING I/O 
I*IOON£ (PUPR) 

IFd.NE.DGo TO 163 

REAO INTO PZ1 *PZ2tPZ3* PRESSURES FOR I- 1 * 64, J=33. 64,Z=1 .64 
I»IRANR <RQ7,PZ1 d»lfl) *98304 ,0 ,PlWR) 

I*IDONE IPUWRI 
IFd.NE.lJGO to 171 

IoIRAnR <RQ7*PZ2(1»1*1) *98304 *192*PLWR) 

I=ID0NE (PLWR) 

IF d » n£*-1 ) GO TO 172 

I*IR«NR (R07,PZ3(l»l»l> *65536 ,384. PLWR) 

wait for outstanding i/o 

I=IDONE (PLWR) 

IFd.NE.lJGO TO 173 
MM»32 
GO TO 36 
CONTINUE 

WRITE PZl»PZ2*PZ3 TO PRESSURES FOR I« 1 *64* J«1 »32*Z«1 * 64 
I*IRANW (RQ7,PZ1(1* 1*1) *98304 .0 ,PlWR> 

I=IDONE (PLWR) 

IFd.NE.lJGO TO 181 

IsIRANW (RQ7.PZ2 ( 1 *1*1) *96304 *192, PLWR) 

I=IOONE(PLWR) 

IFd.NE.DGo TO 182 

I*IRANW (RQ7»PZ3 d *1 *1) *65536 *384*PLWR> 

WAIT'FOR outstanding I/O 
I=ID0NE<PLWR) 

IFd.NE.lJGO TO J 83 

1= I DONE (01 ) ♦iDONEKWj ) *IDONE(W2J*IOONE(W3 ) ♦IdONE(W4J 
IF (I .NE.5J GO TO 191 

#»<**•**«*»#« PHASE III #*»*»**»•*•*•** 

ZH2al 

I0A*1 

I = IOONE1PUPR) * IOONE (PLWR) ♦ I DONE ( DJ ) ♦ I DONE < Wl ) * I qONE {W2J ♦ IDON£ { W3) * 
2 1 DONE ( W4 ) ♦ I DONE ( W8 ) 

IF d*NE»8JGo TO IS'- 

READ PHAT OF PLANES 63 TO 4 INTO PHAT ( 128.64*6) IN LCM 
READ VELOCITIES OF PLANES 63 TO 4 INTO VLG (4?96,3,6 ) IN LCM 
read old velocities of planes 1*2 into vnmi (4096.312) in lcm 
l = lRANR(RQ8.VLGd) *49152, 48iW8) 

IsIDONE ( W8 ) 

IFd.NE.DGO TO 256 
NSECTa4q6 
00 2«5 z=l»6 
I=IOON£(PUPR) 

IFd.NE.DGO TO 203 

I=IDONE (PLWR) 

IFd.NE.lJGO TO 204 

l=I8«NR (R06*PHAT(1,1 ,Z) ,4096*NSECT,PUPR) 

I=IRANR (RQ7 ,PhAT ( 1.2 ,Zj ,4’i96* NSECT.PLWR) 

nsect=nsecT ♦** 

IF (NSECT.GE.512) NSECT=NSECT-512 
CONTINUE 
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5177 DO 40 Z=»l,5 

5201 SMALL IN {PA U ,l),PHAT(i ,1.Z),204S) 

5210 SMALL IN(PA(2049.1) .PHATU .2.Z1.2048) 

5217 SMALL IN(PA{1 ,2) .PHAT (2049,1 «Z) ,2040) 

5226 SMALL IN(PA(2049.2) .PHAT(2n49,2.Z) .2048) 

5235 CALL FFTBXY(Pa) 

C FPTBXY LEAVES THE REAL PRESSURE IN PR(64»64) WHICH IS EQUIVALENCED 
C fO WOLD (1*1.1) 

523b SMALL OUT (P A {1) .PcU.l.Z) .4096) 

5246 4r. CONTINUE 

C LAST DIMENSION FOR PC IS 6 

5230 I0B*1 

5251 DO 50 IPLANEsl *64 

5252 M=* 1 

C TRANSFER ROWS 45 TO 64 TO CALCULATE 47 TO 62 

5253 N j3 A S 

5253 N2=2817 

5254 N3=1280 

5255 N4=l 

5256 NEXT^l 

5257 6n KCZM2-1 

5261 DO 61 Zol,5 

5263 K=k*1 

5264 IF (K.GT .6) K=1 

5267 61 SMALL IN(PD(1»N4«Z)«PC(1»N1*K) «N3) 

5306 SMALL IN (U ( 1 .N4, 3) , VLG (N 6 , l . I0 U ) .N3) 

5320 SMALL IN ( V ( 1 »N4,3 ) . VL6 (N2.2 . 106) .N3) 

5330 SMALL IN<W(1.n4* 3 ) »VLG<N2,3.I0B).n3) 

5340 SO TO (62.63.64.65JNEXT 

5350 62 CALL OVDT (U, V, W. DU*DV«DW.PD*M.CN.CB»cH * RMSU.RMSV.RMSW) 

5367 NU2945 

5370 N2=I 024 

5371 N3=l 

5372 69 SMALL OUT (DU (N3 > . OUlS (N l .1 , IOA ) «N2) 

5404 SMALL OUT (DV (N3 ) *DULG (Nl »2 1 IOA) *N2) 

5413 SMALL OUT (Dw (N3 ) . DUlG (N l .3 . IOA ) *N2 ) 

5422 SO T0(72«73./4,75,76)NEXT 

5433 72 DO 52 Jcl,4 

5435 JJ=J*16 

5436 Oo 52 I si . 64 

5440 00 52 Z«1 *5 

5454 U(I .J.Z)=U<I.JJ.Z) 

5456 V(I.J,Z)=V(I.JJ.Z) 

5^57 W ( I . j.Z)=W(I»JJ.Z) 

5460 52 PO<I»J*Z)=PD<I»JJ*Zl 

C CALCULATE ROWS 63.64 AND 1 TO 14 

5466 N4a5 

5466 Nl*l 

5467 N2= I 

5470 N3=1024 

5471 NEXT*2 

5473 00 TO 60 - 

5473 63 CALL DV0T [U. V«W.DU.DV.DW,PD.M.CN,C8»Cn* RMS U*RMSV.RMSW) 

5512 Nl=3969 

5513 N2=128 
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55H 

5515 

5516 

5517 

5520 

5521 

5522 

5523 

5524 

5525 

5526 

5527 
5530 
5547 

5550 

5551 

5552 

5553 

5554 

5555 

5556 

5557 

5560 

5561 

5600 

5601 

5602 

5603 

5604 

5605 


5605 

5623 

5625 

5627 

5633 

5635 

5636 
5646 
5656 
5673 
5673 
5710 
5710 
5725 
5725 
5742 
5746 

5750 

5751 
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N3nl 

30 TO 69 

73 NEXT“3 
Nl«l 
N2*896 
N3«129 
GO TO £9 

74 Nl*l3 

C CALCULATE ROMS 15 To 30 

N2«769 
N3*1280 
N4*l 

30 TO 60 

64 CALL DVDTfU,V.W.t>U*0V.DW,P0»M»CN,C8»Cll»RMSU.RMsV.RM5W) 

NEXT»4 

Nl»a*7 
N2sl 024 
N3a) 

SO TO 69 

75 Nl=29 

C CALCULATE ROWS 31 TO 46 

N2“1 793 
N3-1280 

N4*>i 

SO TO 60 

65 CALL DV0T(U,VfW.0U»DV.DW,p0«M»CN,C8.Cn»RMSUtRMSV»RMSH> 

NEXT®5 

Nl»l 721 
N2*1024 
N3«l 

GO TO 69 

76 CONTINUE 

C WAIT FOR ANY OUTSTANDING I/O 

r WRITE DUtoVtOW TO- DISK 

c READ PHaT OF IPLANE+4 TO PHATAUOA) 

c REaDVLG OF IPLANE*4 TO VlS(ZM 2 > 

210 loIDONE <W 1 ) *IDONE IW2) * I DONE (W3) *ID0NE <W4) *IOONE (W 8 > 

IF ( I.NE.5) GO TO 21 u 

NSECT=24i»clpLANE4l) 

If insect. GT. 1512) NSECT»NSECT*1536 
IFdTIME.LT. 0)G0 TO 223 
IF ( ISET . tG.O ) GO TO 211 
GO TO (221 *222 *223.224) IM 

211 GO TO (224|22l.22H»223)lM 

22) I=IRANW(RQl«DULS(l»UIOA) ,12288 *NSECT,w1 ) 

GO TO 220 

222 I=IRANW(RQ3*DULG(l*i.I0A) , 12288 . N5ECT.W3) 

GO TO 220 ' 

223 I=IHANW(RQ2,DUL0(l»l.IOA),i2288.NSECT»W2) 

GO TO 220 

224 IaIRANW(RQ4,0ULS{l.i .IOA) , 12288 . NSECT.W4) 

22o IF(M0D(I08»?).NE.0)GO TO ?15 

IF(IPLANE.EQ,64)G0 TO 50 

I0C=I0B-1 

NSEcTaf I PLANE* 4) °24 
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5754 

5760 

5775 

5777 

6003 

6005 

6007 

6011 

6013 

6027 

6044 

6046 

6050 

6053 

6055 

6060 

6070 

6077 

6106 

6115 

6116 
6120 
6123 
6133 
6J34 
6137 


6141 

6167 

6171 

6173 

6174 

6175 
6175 
6205 
6220 
6233 
6246 
6262 
6262 
6275 
6310 
6323 
6337 
6337 
6352 
6365 
6400 
6414 
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IF (NSECT.ST. 1512) NSECT=NSECT- 1536 
l=IRANR(RQ0,Vue(l,l,IOC) , 24576, NSECT,H8) 

215 NSEC 7=8* ( I PLANE* 3) 

IF INSECT. 6T. 504) NSECTsNSECT-512 

217 l=IDONE ( PUPR) 

JFII.NE.DOO TO 217 

218 I=IOONE (PLWR) 

IF ( I «NE. 1 1 GO TO 218 

I=IRANR 1RQ6.PHATAO *1 *104) ,4096*NSECT, PUPR) 

I=IRANR (RQ7 * PHATA (1*2 *IoA) »4096 ,nSECT»PlWR) 

219 I0A=I0A*1 
IF(IPLANE.EQ.64) SO TO 50 
IF ( I0a.gT.2) jOA=1 
108=108*1 

IF ( I0B.GE.5) 108=1 

SMALL IN (PA ( 1,1) .PHATA (1.1. IOA> ,20481 
SMALL IN (PA (2049,1) ,PHAT»(1 ,2.I0a) ,2o40) 

SMALL IN (PA ( 1 *2 ) , PHttTA (2Q49*1 *I0 a) ,2048) 

SMALL IN(PA(2049,2) , PHATA (2049,2, IOA) ,2048) 

C EQUIVALENCE PHaTa(1,1»5) , PHATA ( 1 , 1 , 1 ) 

CALL FFTBXV(PA) 

ZP3=ZM2*5 

IF(ZP3.GT.6)ZP3=ZP3-6 

SMALL OUT (PA ( 1 I ,PC(ltl,ZP3) ,4096) 

ZM2=ZM2*1 

IF ( ZM2.GE, 7) ZM2=] 

5" CONTINUE 
C 
c 

C «*«»*o*ttO***o#« _ PHASE IV ************* 

C 

r 

4(t7 I = IU0NE(Wl),ID0NE(lVi)*I00NE(W3) *IDONE ( W4) + ID0NE ( W0 ) +IDONE (01 ) *ID0N 
? E(PUPR)*IDONE(PLrfR) 

IF ( I ,NE,8) GO TO 407 
IFdTIME.LT.OlsO TO 1000 

if(iset;ne,o>go to 400 

ISETsl 
GO TO 1001 

4q8 GO TO (401*402,403*404) IM 

401 I=IRANR(RQ2, UTNMl (1) ,24576 . 0 ,W2) 

I*IRANR(RQ4,UTN (1), 24576* *W*> 

I=IRANR(RQ3,UN (1) ,24576, 0»w3) 

I=IRANR (ROl *UTNPl (!) *24576, »Wl) 

GO TO 4(j5 

4 U 2 I=IRANR(RQ4,UTNM1 (1) ,24576, 0*W4) 

I = IRANR(RQl,UTN ( 1 ) ,24576, ' ,Wl ) 

I=IRANR(RQ2,QN (1) ,24576, 0,W2) 

I=IRANR(RQ3,UTNPl (1) ,24576, ,W3) 

GO TO 405 

403 I=IRANR(RQl,UTNMl (1) ,24576, 0»W1> 

I=IRANR(RQ3,UTN (1) ,24576, o,W3) 

I=IRANR (RQ4, UN ( 1 ) ,24576. 0 ,W4) 

I=IRANR(RQ2,UTNPI (1) ,24576, *,W2) 

GO TO 405 
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6414 404 I»IRANR{R03,UTNM1 (1) * 24576,0, W3> 

6427 I=IRANR(Rq2*UTN ( 1 ) *24576. «•' *W2 1 

6442 I^IRANR (RQl »UN ( 1 } .24576, 0*Wl ) 

6455 I = 1RAnR(Rq4*UTNP1 (1| *24576* '»**) 

6471 4G5 CONTINUE 

6471 406 I»IDONE(*l! *ID0NE (*2) *ID0NE<*3) ♦ID0NE(W4>*I00NE(W8) 

6507 IF(I'.NE,5)G0 TO 4'0'6 

6511 GO T0{4U,412,413.414)1H 

6521 411 I=IRANR(Rq2*WTNM1 111 *24576*48**2) 

6534 I = IRANR(RQ4,WTN ( 1 )i 24576 *48 * *4 ) 

6547 I=IRANR(RQ3*WN ( 1 ) *24576 *48 **3> 

6562 I=IRANR(RQl.*TNPl (1) *24576*48**1 ) 

6576 GO TO 415 

6576 412 I = IRANR(RQ4.WTNMK1)»24576.48*W4) 

6611 I*lRANR(RQl*WTN 1 1 ) *24576.48*W 1 ) 

6624 Ij=IRANR (Rq2**N ( 1) *245?6*48*W2) 

6637 I = IRANR(RQ3»WTNPltl> *24576*48**31 

6653 GO TO 415 

6653 413 I=IRANR(RQ1*WTNM1 (1) *24576*48**1) 

6666 IeIRANR(RQ3,WTN ( 1 1 » 24576*48. *3) 

6701 I=1RANR(RQ4,WN (1| *24576.48,W4) 

6714 I=IRANR(RQ2,WTNPl (1) .24576.48**2) 

6730 GO TO 415 

6730 414 I=IRANR(R03.WTNM1 (1) *24576*48. *31 

6743 I*IRANR(RQ2.WTN ( 1 ) .24576 *48. *2 ) 

6756 I >* I RANR ( RO 1 * WN (11*24576*48**1) 

6771 I=IRANR(RQ4,*TNP1(1 ) *24576.48.**) 

7005 415 CONTINUE 

7qo 5 DO **5o N=1 * 32 

7007 IF (M0D(N*2) ,EQ, Jl GO TO 430 

7013 IF (N.EQ.l )GO TO 417 

7014 416 l n ID0NE(*lT*t00NE(*2) ♦ I DONE ( *3) ♦ I DONE ( W4) ♦ I DONE ( *8) 

7032 IF(I.NE.S)G0 TO 416 

7034 J=0 

7034 NSECT* (N”2) *48 

7040 GO T0(421. 4*2*423. 424) IH 

7050 421 IsiRANWtRQ2*WN ( I ) *24&76*NSECT«W2 ) 

7064 GO TO 425 

7064 422 I=IRANW(RQ4,WN ( 1 1 *24576*NSECT , *4) 

7100 GO TO 4*5 

7100 423 I3IRANW(RQ1,WN a)*24576*NSECT.Wl) 

7114 50 TO 425 

7114 424 la IRAN* ( RQ3 * WN ( 1 1 *24576* NSECT. W3) 

7l3u GO TO 425 

7130 417 J=1 

7131 425 GO 429 M*=l*6 

7133 Kb (m-1) *4096*1 

7136 Small in(vTnmi (5 i .utnmi (k) .4096) 

7U5 SMALL INtVTN (1),UTN (K)*40g6) 

7153 SMALL IN(VTNPl(l) .UTNPl(K) .4096) 

7161 SMALL IN ( VN (1)*UN (KJ.4096) 

7167 CALL AOVNCtTDT.M.CRMSU,CRMsV,CRMSW,CUMAX.CVMAX*CWMAX) 

7177 SMALL OUT(VNd) .JNlK) .4096) 

7206 428 I=I0QNE<W*)*lD0NE<W2)*IQ0NE ( W3>*rD0NE<W4)*lDQNE(W8> 

7224 IF ( ( (I.NE.5! .AND. (J.EQ.O) ) .ANO* (M,e<3.6) 1 GO TO 428 
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7236 IF< (I.NE.5) .OR, (J.EQ.lHGO TO 429 

7245 J»l 

7245 NSECT«No4B 

7250 SO TO (43l»432»433.434) IM 

7260 43 J IbIRANR(Rq 2,KTNM1 (1) t 24576. NSECT » W2 ) 

7273 I=IRANR(RQ4,lrfTN ( i ) .24576.NSECT. W4 ) 

7306 I«»IRANR(R03»WN ( 1 ) .24576.NSECT * W3) 

7321 I»IRANR(RQ1.WTNP1(1) *24576»NSECT.W1) 

7335 GO TO 427 

7335 432 IslRANR(RQ4,WTNHl (1) .24576. NSECT.W4) 

7350 I=IRANR(R01*WTN ( 1 ) .24576.NSECT.WU 

7363 I=IRANR(RQ2,WN (1) .24576.NSECT.W2) 

7376 I = IRANR(RQ3»WTNP1 (D .24576.NSECT.W3) 

7412 GO TO 427 

7412 433 I = lRANR(RGl,WTNMmi *34576* NSECT.Wj) 

7425 I=IRANR(RQ3.WTN ( 1 ) *24576 »NS£cT .W3 ) 

7440 I“lRANR(RQ4tWN {1) .24576.NSECT.W4) 

7453 13IRANR(RQ2»WTNPH1) *24576*NSECT*W2) 

7467 GO TO 427 

7467 434 I=IRANR(RQ3.WTNMlU)>24576*NSECTtW3) 

7502 I=1RANR(RQ2»WTN ( 1) «24576«NSEcT.W2> 

7515 I=IRANR(RQ1»WN < 1 ) .24576.NSECT.W1) 

7530 I=IRANR(RQ4*WTNP1(1) *24576*NSECT.W4) 

7544 42? CONTINUE 

7544 429 CONTINUE 

7546 GO TO 450 

7547 430 l B IOONE(Wl)+IDONE(W2)*JDONE<W3>*lDONE< w 4>*lDONEtW8> 

7565 IF ( I .NE.5 ) GO TO 43u 

7567 J=0 

7567 NSECTs (N*2) *48 

7573 Go TO ( 441.442.443*444) IM 

7603 441 i*jRANW(RQ27uN ( 1) .2*576. NSECT.W2) 

7617 GO TO 446 

7617 442 I a IRANW (RQ4 • UN (1 ) »24576« N SECT,W4) 

7633 GO TO 446 

7633 443 IsIRANW (RQ1,UN ( 1 ) *245?6*NSECT»Wl ) 

7647 GO TO 446 

7647 444 I*IRANW (RQ3.UN < l ) .24576.NSECT.W3) 

7663 446 CONTINUE 

7663 445 DO 459 M*l*6 

7665 Kb <M*1 ) *4096*1 

7670 SMALL INIVTNHl (1) .WTNM1 <K> *4096) 

7677 SHALL IN t VTN (l).WTN (K ) *4.096> 

7705 SHALL IN < VtNPl ( 1 ) » WTNP1 (K ) . 4096) 

7713 SHALL INIVN (l).WN t K 1 >4096) 

7721 CALL ADVNc<TDTfH.cRHSUfCRMSVtCRMSW,CUMAXiCVMAX»cWHAX) 

7731 SHALL OUT(Vn(1) »»?N<K) ,409 6 ) 

7740 44g I=IDONE<Wl)*lDONE(W2>*ID6NEfW3)*lDONE<W4)*lDONE(W8> 

7756 tF(((I.NE.5). , AND.(J.EQ.On.AND.tM.EQ.6>>GO TO 448 

7770 IF((1.NE.5).OR.(J.EQ.1J>00 TO 459 ' 

7777 J=1 

7777 IF (N.E0.32) GO TO 459 

10002 NSECTaN«48 

10004 GOTO (451 >452,453.454) IM 

10014 451 i=lRANR(RQ2,UTNMl (l) *24576»NSECT,W2) 
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10027 I»IRANR(RQ4,UTN ( 1) ,24576. NSECT. W4) 

10042 Ini PAN R (RQ3,UN ( l ) ,24576 .NSECT ,W3) 

10055 I=IRANR(RQ1,UTNP1 (1) *24576*NSECT»W1 ) 

10071 60 TO 457 

10071 452 I=IRANR(RQ4,UTNMl (1) *24576. NSECT, W4) 

10 1 04 I*IRANR(RQl,UTN ( 1 ) .24576*NSEcT.Wl ) 

10117 I«IRANR(RQ2,0n fl) *24576. NSECT. W2) 

10132 I*=IRANRIRQ3»UTNP1 (1) ,24S76*NSEcT.W3) 

10146 GO TO 457 

10146 453 IcIRANR|RQ1*UTNH1 (1) *24576. NSECT, WH 

10161 I=IRANR(RQ3,UTN ( 1 ) ,24576*NSECT»W3) 

10174 I=IRANR(RQ4,UN (II *24376. NSECT, W4) 

10207 I=1RANR(Rq2,0TNP1 (1| ,24576*NSECT,W2) 

10223 GO TO 457 

10223 454 l=IRANR(RQ3»UTNMl (1) ,24576, NSECT. W3I 

10236 1=IRANR(RQ2*UTN ( ll *24576*NSECT,W2) 

10251 I*IRANRfRQl t UN (1 ) ,24576*NSECT.W1 ) 

10264 I=IRANR<RQ4,UTNPI (1) , 24576, NSECT. W4) 

10300 457 CONTINUE 

10300 459 CONTINUE 

10302 450 CONTINUE 

1 0304 GO T0(5ll, 512*513. 514) IM 

10314 511 IriRANW (RQ2,WN { 1 ), 24576* 1488 . W2) 

10330 GO TO 515 

10330 512 I*IRANK(RQ4,WN ( 1 ) *24576 . 1488, K4) 

10344 G 0 TO 515 

10344 513 I=IRANW(RQ1,WN ! 1 ) ,24576, 1488. Wl > 

10360 GO TO 515 

10360 514 I = IRANW(RQ3,WN 1 1 1 *24576* 1488.W3) 

10374 515 CONTINUE 

10374 1000 ITIM£»ITIME*1 

10376 IF (I tlHE.EO.O ) GO TO 1001 

10377 ISeT^O 

10377 RMSDsSQRT (RHSO) /512, 

10403 RMSU=SQRT(RMSU)/512. 

10407 RMSV=SQRT(RhSV)/512. 

10413 RMSW»SQRT(RMSW)/S12. 

10417 UMAX=SQRT((JhfAX) 

10421 VMAX»SQRT(VMAX) 

10423 WHAX“SQRT(WMAX) 

10425 CRMSU*SQRT(CRMSUI/S12, 

10431 CRMSV=SQRT(CRHSV) /512, 

10435 CRhSW=SQRT(CRMSW) / 5l2. 

10441 CUMAXsSQRT(CUHAX) 

10443 CVHAX=SQRT(CVMAX) 

10445 CWMAX=SQRT <CWMAX) 

10447 SKlaSKl/SQRT!SK2<H>3)*12B, 

10456 PRINT 899,RhS0*R4SU*RhSV*rmSW*Ski 

899 FORMAT!* RMSD»*E‘l 4, 7* RMsUb*E 14,7* RMSV=*El4.7* RMSWn*El4.7* 
2SKEWNESS= a *El4.7) 

10473 PRINT 880*OmAX,Vv)AX,WMAX 

888 FORMAT!* UMAX**Eil4.7* VMAX=*eU, 7 * WMAXn*El4.7) 

10505 PRINT 895 » CRMSU* CRMS V * CRMSW 

10517 PRINT 888*CUMaX*CVMaX*CWMAX 

895 FORMAT (* * RMSUa«El4. 7* RMSV"*E14.7* RMSW**E14,7) 
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10531 

10533 

105*1 

105*2 

10543 

10546 

10552 

10560 

10562 

10570 
1 0600 
10600 

10604 
1 0605 
10607 
10623 
10625 
10627 
10631 
10634 
10643 
10650 
10654 
10657 
10661 
10675 
10675 

10701 

10702 
10704 
10720 
10 722 
10724 
10726 
10731 

10740 
10745 
10751 
10 754 
10756 
10772 

10775 

10776 
11001 
11002 

11005 

11006 
Hull 
11012 
11015 
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XXa20./(3.l4l5927»64.«15360. ) 

00 889 Icl,64 
SPI(I)=0. 

889 SPR(I)*XX*SPR(I) 

CALL F 7T2 (SPR <1 1 tSPl (1 ) *64»-l ) 

PRINT 887 

887 FORMAT {* Ell=a) 

Print 886*ispru> «i»i*32) 

886 FORMAT [ 1X*8£12,5) 

CAll"SECOND'(TX) 

PRINT 1900 »TX 

1900 FORMAT (p CpU TIME =*E14,7) 

GO T0(l500ii200»l60o*130j) iM 

1200 CONTINUE 

IF t IDONE (W4) ,NE. l ) GO TO \z 0 0 
C SAVE'TIME STtp 3 FROM W4 IN FSET7 
DO 1210 K=1 * 8 
NSECT*(K"1)°192 

IaIRANR(RQ3*DUMl (11 *98304, NSECT*W4) 

1201 I=ID0NE(W4) 

IF(I.NE,1)G0 TO 1201 
00 1210 J=l»3 
JK=(J-1)«32768+1 

SMALL IN ( DUMSM ( 1 1 *DuMl (JK) ,32768) 

WRITE(7) (OUMSM(I) ,1=1*32768) 

1210 CONTINUE 

R=(ITIME-40) /2 
END FILE 7 

GO TOI1401 *1 402* 1403 *1404, 1405 *1406*1 Ao7* 1408) M 

1300 CONTINUE 

IF ( IOONE ( W3) jNE, 1 ) GO TO 1300 
C SAVE TIME STtp 5 FROM W 3 IN FSET 7 
DO 1310 K=l,i 
NSECT=(K-U»l92 

I*IRANR(RQ3«DUMl (1) , 98304, NSECT,W3) 

1301 i=io6ne<w3) 
iFd.Ne.nao to 1301 
DO 1310 Jsl *3 

JK= < J-l ) «32768*1 

SMALL IN ( OUMSM ( 1 ) »QUMl ( JK) ,32766) 

SRlTE (7) (oUMSM(I) *I»1*327&8) 

1310 CONTINUE 

M«=(ITIME-40)/2 
END FILE 7 

GO TO (1401 *1402* 1403 * 1404, 1405, 14o6*1407* 1408 )M 

1401 CALL AFSR£L(5LFSET7*0,I1TAPE) 

GO T01269 

1402 CALL AFSREL(5LFSET7*0*I2TAP£) 

So TO 1269 

1403 CALL AFSREL(5LFSET7,0*I3TAPE) 

So TO 1269 

1404 CALL AF$REL(5LF5ET7»0»I4TAPE) 

GO TO 1269 

14 0 5 CALL AFSREL(5LFSET7*0*I5TAPE) 

GO TO 1269 
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11016 

11021 

11022 

11025 

11026 

11031 

11032 
11032 
11036 
1104-6 
11050 
11052 
11066 
1 1070 
11072 
1 1 0T4 
11077 
11106 
11113 
11117 
11120 
11120 
11124 
11134 
11136 
11140 
11154 
11156 
11160 
11162 
11165 
11174 
11201 
11205 
11210 

11236 

11240 

11250 

11252 

11254 

11270 

11272 

11274 

11276 

11301 

11310 

11315 

11321 

11323 

11325 

11341 

11343 

11345 

11347 
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1406 CALL AFsREL(SLFSET7,0.I6TAPE) 

80 TO 1269 

1407 CALL AFsreL(5LFSET7,0iI7TAPE) 

GO TO 1269 

1408 CALL AFSREL(5LFSET7,Q,I8TAPE) 

GO TO 1269 

1500 CONTINUE 

IF(IDON£(W2| ,NE.l)GO TO 1500 

CALL OPEN (5lPSET7» 0.2340 006 *0,0. 0*3200. IT APE) 

DO 1510 K»l,a 
NSECT*(K-l)*l92 

InJRANR (RQ2.DUM1 (U .98304. NSECT .W2) 

1501 I*ID0NE (W2) 

IF ( I .N£» 1) GO TO 1501 
DO 1510 J“l»3 
JK« 1 J"1 ) *32768* 1 

Shall inidumsm(I) .duhkjkj .32768) 

WRITE (7) <OUHSM(I) »I«l*327i8) 

1510 CONTINUE 

GO TO 1269 
1600 CONTINUE 

IF (I DONE (Ml I .NE.l)GO TO 1600 

CALL OPEN {5lESET 7 » 0. 23400 00 .0*0*0. 3200. ITAPE) 

DO 1610 Ko l.s 
NSECT=(K-l)*l92 

I=IRANR(RQ1*DUH1 (1) 198304*NSECT*W1) 

16C1 I*I0ONe(W1) 

IF 1 1. NE.DGO TO 1601 
00 1610 J=r,3 
JK=(J-1) 032766*1 

SMALL IN (OUHSfttl ) .DLJMl (JK) .32768) 

WR JTE ( 7 ) (PUMSM(I), Ini, 32769) 

1610 CONTINUE 

1269 IF(ITIME.LT»52)G0 TO 10O1 

1100 I a IDONE (Dl)*IDONE(WD*IOONE(W2)*lDONEtW3>*lDONE(W4)*IOONE(W8) 

2*ID0NE(PUPR) ♦lOONE(PLWR) 

IF 1 1 ,NE,8) GO TO 11 Ou 

CALL OPEN (5LFSET5.0 .2340008* 0» 0*0 .3200 .ITAPE) 

DO 1120 ' KbI .8 
NSECT® (K-l ) *192 

InIRANR (RQltOUMl(l) , 98304, NSECT » W1 ) 

1102 I=IOONE(Wl> 

IF(I. NE.DGO TO 1102 
00 1120 Jnl ,3 
JK«(J-1 ) *32768*1 

SMALL INtDUMSM(l) ,DUMl < JK) .32766) 

WRlTEtS) (OUMSM(I) ,Jal. 3270s) 

H2l> CONTINUE 

00 1110 K-l, 8 
NSECT* (K-l 1*192 

l5=IRANR(RQ2»DUMl (1) .98304.NSECT.W2) 

1101 I^IOONE <W2) 

IF Cl • NE.DGO TO 1101 

00 mo jnl . 3 
JK=(J*1 1*32768*1 


104 



RUN-LCM89 GAP 74/11/14 16.32.37 VSClARIUH* pA6t NO. 22 

11352 SMALL IN ( DUMSM ( 1 ) ,OuMl (JK) ,32768 1 

11361 WRITE (5) [DUMSM(I) ,1*1,32768) 

11366 llln CONTINUE 
11372 END FILE 5 

11374 CALL AFSREL(5LFSET5»0»I05TAPE! 

11377 CALL OPEN15LFSET6.0. 23400,9.0*0*0. 3200*ITAPE) 

11407 DO 1140 K=l,8 

11411 nSECTc (K' l I *192 

11413 I=IRANR(RQ3.DUM1 (1) .98304, NSECT.W3) 

11427 1104 I=I00NE(W3) 

11431 lFtI.NE.U60 TO 1104 

11433 DO 1140 Jsl . 3 

11435 JK*(J-U*32768.1 

11440 SMALL INtDUMSMtU.OUMltJK). 32768) 

11447 WRITE (6) (oUMSM(I) ,1*1*32768) 

11454 1140 CONTINUE 

11460 60 1130 Kcl.8 

11462 NSecTs (K" l ) °192 

11464 I*IRANR(RQ4,DUMHU *98304, NSECT*W4) 

11500 1103 I*=IoONE'tW4) 

11502 IF(X.NE,1)G0 TO 1103 

11504 00 1130 J*l,3 

11506 JKb(J-I) *32768*1 

11511 SMALL iNtDUMSM(l) .OUMKJK) ,32768) 

11520 WRITE (6) (DUmSM(I) .1*1,32768} 

11525 1130 CONTINUE 

11531 END FILE 6 

11533 CALL AFSREL (5LFSET6 * 0. 1 06TAPE) 

11536 00 1131 1*1,5 

11544 I05TAP£(D*I07TAPE(I> 

11545 1131 I06TAP£(I)*I08TAP£(I) 

11547 IF(ITIMe.eQ.24)GO TO 1001 

11551 STOP 

11553 END 


PROGRAM LENGTH INCLUDING I/O REQUEST TABLES - GAP 
15131 


statement 
STMT NO* 

ASSIGNMENTS 

location 

stmt 

NO* 

location 

STMT 

NO* 

location 

STMT NC 

13 

* 

3563 

16 

* 

4057 

17 

ft 


19 

20 

f* 

4302 

22 

* 

3651 

23 

r* 

3726 

24 

25 

l* 

3774 

26 

ft 

3547 

28 

ft 

3556 

36 

38 

* 

4760 

5 J 

rt 

6140 

60 

ft 

5260 

61 

62 

r* 

5351 

63 

r* 

5474 

64 

ft 

5531 

65 

69 

* 

5373 

72 

1* 

5434 

73 

it 

5517 

74 

75 

* 

5555 

76 


5606 

100 

rt 

2406 

109 

113 

* 

4122 

116 

* 

4221 

117 

f* 

4225 

118 

130 

ft 

1173 

131 

i* 

1220 

132 

it 

1247 

133 

134 

* 

1 3?5 

139 

r* 

1354 

149 

ft 

4304 

150 

IS] 

r* 

1441 

152 

T* 

1455 

153 

it 

14?1 

154 

155 

* 

1521 

156 

ft 

1523 

157 

it 

1540 

161 

162 

<t 

4656 

163 

ft 

4675 

171 

(♦ 

4714 

172 
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SUBROUTINE flOVNC ( TDT.M *RMSU* RMS V.R hSW. UMAX* VMAX.WMAX) 


20 


COMMON U ( 32768) 


C 

CU1./(2.*DELT> 

20 


Cl=i./TDT 


C 

C2=1./(2.*DELT**2) 

21 


C2»2./T0T*«2 


C 

C3=DELT 

23 


C3».5*TDT 


C 

C* b 0ELT#*2/H, 

25 


C*b. 125*TDT»«2 v 


c 

C5*DELT*»3/6, 

27 


C5=T0T**3/2*. 

30 


C6»C**Cl 

32 


C7n C 5«C2 

35 


DO 10 I»1 * *096 

A3 


U ( I) =U(H*C3*U (1*8192 > *C6* (Ut I* 12288) »U( 1**0961 > *C7* (U < I* 12288) - 



22.«U(I*8l92 ) *U( I **096) ) 

53 


10 CONTINUE 

5* 


GO TO(2o»3o»40«20*30**0)M 

71 


20 DO 28 JbI * * U 96 

73 


X=U(I)«*2 

75 


Rmsu=rmsu*x 

76 


If (x.gt.umaxjumakbx 

101 


25 Continue 

103 


GO TO 5o 

1 0 A 


3n DO 38 l3l»*096 

106 


Xau( I | **2 

uo 


RMSV»RMSV*X 

m 


IF(X.GT.VMAX) VMAXsX 

u* 


35 Continue 

116 


So TO 5 o 

117 


*0 DO *8 I»lt*o96 

121 


X=(J(II »*2 

123 


SMSW»RMSW*X 

12* 


If (x.gt.wmaxiwmaxsX 

127 


*s Continue 

131 


5» RETURN 

132 


ENp 


SUBPROGRAM LENGTH - ADVNC 
162 

STATEMENT ASSIGNMENTS 


stmt NO* LOCATION 

STMT 

NO* 

LOCATION 

stmt 

NO* 

location 

STMT Nf 

20 * 72 

25 

* 

102 

30 

* 

105 

35 

*0 * 120 

BLOCK NAMES AND LENGTHS 
•» 100000 

*5 

* 

130 

50 

* 

132 


VARIABLE ASSIGNMENTS 

name * LOCATION 

NAME 

* 

location 

name 

* 

LOCATION 

NAME 
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SUBROUTINE DVRGNC(C12,RMS0.«MSU»RMSV,RMSW»SK1.SK2,UMAX,VMAX.WMAX,M 

?« ID.OEUT ) 

25 L.CM/8B1/VLG 

25 LCM/882/WLG 

25 LCM/BB3/0U6 

25 DIMENSION VUG (4096*3. 6} »WLG(4096*3,8) .DUG ( 64 , 64 . 16) .A(20 480). 0(819 

22) 

25 COMMON DUMSM (32768* 

25 EQUIVAUENCE t A ( 1) . DUMSM ( 1 ) ) * (D < 1 ) . DUMSM < 20481 ) ) 

C DISPLAY OFP 

25 N»M0D(M,4) 

31 IF (N.EQ, 0) N=4 

33 0T0T=>D£UT/12. 

35 TTDT»2. o 0EUT/3. 

40 DO 5. K=3,6 

41 IF{(N.EQ.2).0R.(N.EQ.^))G0 To 2 

50 SMALL IN(A(3) ,VLG(1 .l.K) ,4096) 

57 GO TO 3 

57 ? SMALL IN(A(3> *WLG(1.1*K) ,4o96) 

67 3 A ( 1 ) =A (4097) 

71 A (2) =A (4098) 

72 A(4p99)aA(3) 

74 A ( 4 1 00 ) =A ( 4 ) 

75 00 10 1*1.4096 

104 10 D ( I ) 0 A(I)-A(I*4)+8.»(Ad*3)-A(l+l)) 

111 IF ( (N.EQ.2) .OR, (V^EQ.4) )GO TO 12 

124 SMALL IN(A (129) » VL 8 (1.2?K),4096) 

133 GO TO 14 

133 13 SMALL IN ( A ( 129) . iJLG ( 1 .2 .K) . 4o96> 

143 14 00 IP 1=1,128 

150 15 A ( I ) °A ( I *4096 ) 

152 DO 16 1=4225.4352 

162 16 A(I)=A(I-409S) 

164 DO 2*' 1=1,4096 

176 20 D ( I 1 *0 ( I ) ♦ A(I)«A(I*256)*8.<MA(I*192)-A(I*64>> 

204 J=K- tf 

205 JJ*1 

206 IF ( (N.EQ.2I .OR. (M.E0.4) ) GO To 23 

222 00 26 1=1,5 

223 SMALL IN ( A < JJ> , VLG < 1 »3. J) ,4096) 

233 J=j*l 

234 25 JJ=JJ*4096 

237 GO TO 24 

237 23 DO 2t> 1 = 1.5 

241 SMALL IN(A(JJ) ,WtG<) »3»J) .4096) 

251 J=J*1 

252 26 JJ=jJ*4096 

255 24 00 30 1=1,4096 

263 30 D(I)»ClZo(D(I)*A(I)»A(I*16384)*B.»(A(I,122B8)-A(I*4096))) 

271 J=K-** IN-1 ) "4 

275 GO T0(40,32,33) ID 

310 32 SMALL IN(D(4097) ,DLG(1,1,J) ,4096) 

320 00 36 1=1,4096 

325 Q(I)*'»0TDT 6 D(I)*D(I*4096) 

327 36 CONTINUE 
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QVRGNC 


74/1 1/1* 


16,32.37 V5CLARK1HA PAGE NO. 2 


330 50 TO 4o 

334 33 SHALL IN(D(4097).OLG<1.1,J)*4096> 

344 gO 37 jol ,4096 

351 0tI) B TTDT*0(I) *0(1*4096) 

353 37 CONTINUE 

357 4n SMALL OUT <0 ( I > .DIG < l * 1 tJ> * 4096) 

367 50 CONTINUE 

371 RETURN 

372 END 


SUBPROGRAM LENGTH - DVRGNC 
424 


statement 

5TMT NO* 

assignments 

LOCATION 

STMT 

NO * 

location 

STMT 

NO* 

LOCATION 

STMT NO 

2 

* 

60 

3 

& 

70 

13 

r* 

134 

14 

23 

r* 

240 

24 

& 

256 

32 

r* 

311 

33 

40 

& 

366 









block names and lengths 

r» lOOOOO 


VARIABLE ASSIGNMENTS 


NAME 

i* 

LOCATION 

NAME 

i* 

location 

NAME 

f* 

LOCATION 

NAME 

A 

i*R 

OCCH 

CI2 

MR 

0 

D 

r*R 

50000C01 

DELt 

UUMSM 

*R 

0C01 

I 

(•I 

415 

10 

,*1 

13 

J 

JJ 

*1 

41 7 

K 

r»I 

420 

M 

r* 1 

12 

N 

OTDT 

*R 

422 

rmsd 

i*R 

1 

RMSU 

r*R 

2 

RMSV 

KMSN 

*R 

4 

S<1 

i»R 

5 

Sl<2 

«*R 

6 

ttdt 

UMAX 

*R 

7 

VMAX 

i*R 

10 

wmax 

r»R 

11 


LCM BLOCK 
BB} * 

NAMES AND LENGTHS 
300000 

682 

i» 

300000 

BB3 

4 

200000 


LCM VARIABLE ASSIGNMENTS 
Name * LOCATION 

Name 

i* 

location 

name 

4 

location 


ULG 

*R 

0L03 

VLB 

i»R 

OL01 

WLG 

r*R 

0L02 



external assignments 

ACGOER r*R 


START OF - CONSTANTS TEMPORARIES INDIrECts • UNusEO COMP I 

374 403 413 733 
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"UN-COMBO 


74/11/1* 


16. 32. 37 


V5CLAHMHA PAGE NO. 1 


SUBROUTINE FFTFXY[P) 

3 DIMENSION P(A 0 96*2t 

3 DO 1 I*lt*0g6 

7 1 P | I , i ) =0 . 

11 CALL FFT2fPa*l) tpliiZI .*096,-1) 

15 RETURN 

16 END 


SUBPROGRAM LENGTH - FFTFXY 
30 

VARIABLE ASSIGNMENTS 


name * 

LOCATION 

NAME 

r* location 




i *1 

27 

P 

<*R 0 




EXTERNAL 
F FT2 *R 

STAPT OF 

assignments 

CONSTANTS 

20 

temporaries 

22 

INOIrECTS 

24 

- 

UNUSED COMPI 
7*6 
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74/11/14 


16.32.37 


V5CLARK1HA PAGE N 0 . 1 


SUBROUTINE FFTBXY (P) 

3 DIMENSION P(4096*2) 

3 CALL FFT2(P(Itl) »PU*2) i4096»1) 

7 Return 

10 END 

CD Fus 


SUBPROGRAM LENGTH - FFT8XY 
20 

VARIABLE ASSIGNMENTS 
NAME •* LOCATION 

P i*R 0 

EXTERNAL ASSIGNMENTS 


F FTg *R 

START OF 


constants 

TEMPORARIES 

INDIrECTS 

• 

UNUSED COMPI 



12 

14 

16 


746 
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16.32.37 V5CIARK1HA 


PAGE NO. 1 


SUBROUTINE DELSQ(P»CtC64,S.S64*l»N,C6> 

20 DIMENSION P(256,S4,2) ,C(25Si ,S(255j ,C64 (255) ,564(255) 

C DISPLAY OFF 

20 DO 10 Msi,256 

21 CALL FFTH(P(M*1»1) iP (M*lt2) *64»-256) 

34 in Continue 

36 00 20 Ksi.64 

37 K2»2*K-1 

40 K3»3»K-2 

42 K4=4»K-3 

45 UO 2V Jsl *4 

46 Jl*j+L*N-1 


50 J2»2*J1-1 

53 j3a3*>jl-2 

55 J4=4« Jl«3 

57 DO 20 Jwl»64 

61 |2=2»I-1 

62 l3n3*I-2 

64 I4a4»l»3 

66 X=C6»(2. <, (C(I4)OC64(J4)-S(I4)«S64(J4) *C64 i 14) *C64 <K4> ) *128.* (C< 12) 

2*C64 { J2 ) -sl l2 ) *S64 ( j2 ) *C6? < 12 > ♦ C64 («2 ) ) -32 . * { C ( 13 ) *C64 { j3 ) -S ( 13 ) * 
3S64 (J3)*C64 (I3J+C64 (K3) ) *32.*(f< I) «C64 ( Jl ) -S ( I ) *S&4 ( Jl ) *C64(I) ♦ 
4C64<K) )-390.) 


151 



f J»l)#64 

155 



IFUPS(X) .LT.1.E-05JGO TO 17 

161 



X*l./X 

162 



GO TO 18 

163 


17 

X=«0. 

171 


18 

P(M,K,1)*x«P<M,K,1) 


C 


DISPLAY ON IFUK.EQ.ll.AND.U. EQ.l)) 

172 



P(M.K.2)aX*p(M.K.2> 


C 


DISPLAY off 


C 


IF t (K.EQ.l) .AND. ( I.ES. 1) ) PRINT 51 »K* I *P (Mt K. 1 ) 


C 

51 

FORMATS P (M,K » 1 ) *«*El4« f ) 

176 


20 

CONTINUE 

204 



DO 30 M*l» 256 

206 



CALL'FFT 2 (P(M.l.l) .P(M.1,2) .64,256) 

221 


30 

CONTINUE 

223 



RETURN 

224 



END 


SUBPROGRAM LENGTH - DElSQ 
300 

STATEMENT ASSIGNMENTS 


STMT NO* LOCATION 

STMT 

NOi» 

LOCATION 




17 * 164 

18 

* 

165 




VARIABLE ASSIGNMENTS 
NAME * LOCATION 

name 

<* 

LOCATION 

NAME 

* 

location 

L r*R 1 

C6 

*R 

7 

C64 

*R 

2 

12 *1 262 

13 

*1 

263 

14 

*1 

264 


NAME 

I 

J 
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SUBROUTINE UVDT ( J'» V .W« Du, DV * OW » P.M, CN» CB, Cll .RMSU.RMSV. RMSH) 

26 DIMENSION U(64oo> *V(6400) *W (64001 »DU ( 1024) * DV ( lo24> ,DW(1Q24) t 

2p (6400 ) 

26 00 10 I a 1 . 1 024 

56 DU (I)«U( 1*2688) -C8* (P< I*26B6!*P<I *2690) *8.* (P(I*26B9)-P( 1+2687) ) > 

66 DV (I ) sV ( I +26BB) -C8* (P ( 1+2560) ”P( I +2816) +8.*rP(I+ 2752) -P( 1+2624) ) ) 

77 10 DW< I ) «=W(I + 268B)-C6*(P< 1 + 128) -P (1+5248) +8,* ( P ( I+3 96B) *P( 1 + 1408) ) ) 

141 RETURN 

142 END 


SUBPROGRAM LENGTH - DVDT 
203 

VARIABLE ASSIGNMENTS 


NAME 

4 

LOCATION 

NAME 

i+ 

LOCATION 

NAME 

f* LOCATION 

Name * 

ON 

*R 

10 

CU 

r+R 

12 

C8 

r+R 

11 

Du ■) 

L)V 

j+R 

4 

Dri 

r+R 

S 

I 

«+I 

202 

M + 

P 

+ R 

6 

RMSU 

r+R 

13 

RMSV 

i+R 

14 

RMStf 

u 

+ R 

0 

V 

i+R 

1 

W 

-*R 

2 


start 

OF 

- 

constants 

144 


TEMPORARIES 

145 

INDIRECTS 

163 

- 

UNUSEO COMP 11 
743( 


112 
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SUBROUTINE CAUCPR(U,V.W*P,M»C3,C7) 


17 


OIMENSION U (6400) ,V(6400) ,W(6400) ,p(4096) 

C 


DISPLAY OFF 

17 


If (M.EQ.2) 60 TO 20 

20 


11 = 1 

21 


12=1024 

22 


IF(M.EQ.I) J=2944 

25 


IF(M.EQ.3) JP896 

30 


1F(M,EQ.4) j=l92r 

33 

5 

00 10 l jb 1 1 , j2 

37 


JJ=I»J 

40 


P < JJ) sC7*P < JJ> *C3* (U ( I*2686> -U (I*269o)*V( 1*2560 )-V( 1*281 6) *W(I*1 28 

?) -H (1*5248) *8,0 1+5687) ♦ V ( I 752) - V (1*2624) +W (1*3968) 

3-W( 1+1408} ) ) ' 

65 

Id 

CONTINUE 

67 


IF (M.NE.2) GO TO 50 

74 


GO TO 3„ 

74 

2 1 

1 1 = 1 

75 


12=128 

76 


J*3968 

77 


GO TO 5 

lOo 

30 

IF(J.EQ. - 128)60 TO So 

102 


11 = 129 

103 


12=1024 

104 


J=~l 28 

105 


GO TO 5 

106 

So 

RETURN 

107 


END 


SUBPROGRAM LENGTH - CALCPR 
142 

STATEMENT ASSIGNMENTS 


STMT 

NO* LOCATION 

STMT 

NOf* 

location 

STMT 

NOi* LOCATION 

STMT NO* 

5 

1* 34 

20 

l» 

' 75 

30 

<* lol 

50 * 

VARIABLE ASSIGNMENTS 







name 

* LOCATION 

name 

I* 

LOCATION 

NAME 

* location 

name * 

C3 

i»R 5 

C7 

i*R 

6 

I 

+ 1 135 

n * 

12 

*1 137 

J 

.+1 

140 

JJ 

,•1 141 

M 4 

P 

+ R 3 

u 

<*R 

0 

V 

,R 1 

w + 

start 

OF 

CONSTANTS 


temporaries 

indirects 

unused comp it 



111 


112 


120 

743i 


113 



37 

37 

41 

43 

46 

67 

107 


224 


356 


537 

541 

541 

563 

603 


RUN-LCM89 74/11/14 16,32.37 v5ClARK 1HA pA6§ NO. 1 

SUBROUTINE DVDTMPi(U,V*W.DUiOV.DW.M,CN.CeiCll.C12.ITlME.RHSD*RMSU» 

2RHS V.RMSW, SKI, SK2, UMAX, VMAx.HHAx.S.IsET) 

DIMENSION U (64,80,5) ,V (640o ! • W^OO* ,OU( 1024) ,DV < 1024) ,DW(1024> , 

2P (6400) »S (64) 

Cl=.5 tf CB 

C2»Cll*CN 

I F t ISET.EQ, 1) 50 TO 20 
UO 1« I»1 » 1 024 

DIV=U( 1*2686 ) -U( 1*2690) +V( 1+2560)- v(I*28lb) *W( 1*1 28)-W( I* 5248) 

9*8. » (U ( I *2689) -u ( I *2687) *V ( 1*2752) -V( 1*2624 )*W< I *3968) -W( 1*1408) ) 

UU( I )3-Cl*(U( 1*2688) *(U ( 1*2686) -U 1 1 *2690) *8. * ( U ( I *2689) -U { 1*2687) ) 

2 f + V ( I + 2688)o(U(I*H560)-U(V*28l6) *8 . * ( U ( I +2752 ) -U 1 1 *2624 ) ) j 
3*H ( 1*2688) * (U ( 1*128) -U (1 + 524 8) *8.* tU ( 1*3968) -U < 1*1408) I ) 

4*U( I *2686) *o2-ud*2690) »»2 + 8,*((J ( 1*2689) 6*2-U( I +2687) ***2) 

5*1) ( 1 + 2560) *V ( I*2560)-U(I*28l6)*V (1*28) 6) *8,» (U < 1*2752) *V (1*2752) 

6-Uf 1*2624) *V (1*2624) ) *U t j*1 2$) *W ( 1*128) -U ( 1*5248) *W (1*6248) *8.* 

7(U( I+3968)<>H(I+3968)-U(I* i 408)*W{I*l408) )*U(2688)*DIV) 

9*C2* ( 16, * (U ( J *2689) *U ( 1*4687) *U < 1 *2752 )+U( 1*2624) *U( I +3968) *Ut I* 

6 1^08) ) -0(1*2686) *U( 1*2690 )-U( I *f 81 6) -U ( 1*2560 ) ~u ( 1*128) -U (1*5248) 

B-90 , a U ( 1+2688) ) 

UV( I) »-Cl«7U (1+2688) «tV(l*2686)-V(I*269O)+0.»(V(I + 2689)-V( 1*26871 ) 
?i*V(I+2688!«(V(I*2560)-v (1+2816) *8.* ( V ( I *2752 ) -V ( I *2624 ) > ) 

3*W( 1+2688) * (V( I* 12B)-V( 1*5248) +8.* (V ( 1*3968) -V 11*1408) ) ) 

4*V ( 1+2686) «U (I *2686) -V ( 1+2690) *U (1*2690) *8.* (V ( 1*2689) *U (1*2689) 

5- V( 1*2687) »Ut 1*2687) ) *V < 1*2560 ) **?-V < I*2«l 6) **2*8.* (V <I+?7S2> «»2 

6- V(I»2624)»»4)*V(I*if8>«W<]: + 12$!-V(I*5248) <*W ( 1*5248) *8. <MV (1*3968) 
7*W(I*3968)-V(I*140B)*W(I**4 8) ) + V (2688 > *DI V ) 

9+C2* (16.4 (V ( 1*2689) *V(I+ C 687)*V( 1+2752) *V( 1*2624) *V( 1+3968) +V( I* 

A1408) )-V( 1*2686 )-VU*2690)-V (I*f81 6) -V( 1*2560 )-V( 1*128) -V( 1*5248) 

B-90 .»V(I*2688) ) 

OW ( I) =«cl °I(J ( 1*26881* (W ( 1+2686) -W( 1*2690) *8. * (M ( 1*2689) -W ( I + 26B7) ) 
2)+V(I*2688)*(W(I*2565 !)-W(t+ 2816) *8, * ( W < I *2752 ) -W ( I *2624) j ) 

3+W( 1*2688) *(«(!+ 128) -W (I*=?4b) *8.»(W( 1+3*68) -W( 1*1408) ) ) 

4* W ( 1*2686) *U ( 1*2686) -W ( 1*2690 )*U (1*2690) +8, *(W( I *2689) *U(I*2689) 

5- W ( 1*2687) »U ( 1*2687 ) ) ♦ W < 1*2560 | «v ( 1+2560 ) ( 1*281 6) «v( 1+2816) +8,* 

6 ( W C 1*2752) *V< 1*2 752 )-W( I *2624 1 *V (1+2624) ) ♦ M ( J*1 28) **2-W ( 1*5248) **2 
7+8.*(W(I*3968)**2~W(I+j:408)**2)*W(2688)*UIV) 

9+C2* (1 6. * (W ( 1*2689) *tf ( I *2687) *W( 1*2752 )*W( 1*2624) *W( 1*3968) +W( I* 

4*40 8) ) -W ( I *2686) -W ( 1*2690 ) -H ( l*fa) 6 1 -W <1+2 560 )-W( 1*1 28) -W( 1*5248) 

H-90 . *W ( I *2688 ) ! 
lt> CONTINUE 
80 TO 5p 

20 DO 23 Ul.1024 

DI V a U (1*2686) -U( 1+2690 ) *y < 1*2560 ).v(i*20l6)*tf(l+128)-W (1*5248) 

2*8.* (U ( 1 *2689) -U ( I *2687) *V ( l *2752) -V( I *2624) *W( 1*3968) -W( I *1408) ) 

DU( I) =-Cl* (U U *2688) * (Ut 1*2686) *U ( 1 + 2690) + 8, * (U ( I *2689) -U ( 1*2687) I 
2) *V (1*2688) *(U(I+2560)-U( 1*2816) +0.*(U( 1*2752) -Ut 1*2624) ) ) 

3*W ( 1 + 2688)* (U( I* (28) -U( I + 5?48) *8.» (U ( 1 + 3968) -U (1*1408) ) ) 

4+U ( 1*2686) **2-U 1 1+2690) **2+8, *(U (1*2689) **2-U ( 1+2687) **2) 

5*U( 1*25601 *V( 1*2560) -U( I *2816) »V ( i*2{jl6) +8. » ( U ( I + 2752) *V( 1*2752) 

6- U (1+2624) *V (1*2624) ) *U (I* 1 28) *rf ( I*l28)-U ( I *5248) «W (1*3248) +B.» 

7(U( J*3968) «W< 1*3968) -UC J* 1408) »W ( J* 14 08) > *(J (2688) *OlV) 

9*C2*(l6.*(U(I+2b89) +U ( 1 +2 6 87) *U 1 1 + 2752) +U (I *2624 )+U (1*3968 1 *UU* 

Al4ns) ) -U ( 1+2686) -U ( I +269Q) -U ( I*f8l6) -U( 1*2560 )-U( I* 128) -U( 1+5248) 
8-90.4U(I*2688) ) 


REPRODUCIBILIXV Of ‘ 
0RIGMAL'PAOB II F90& 
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720 0V(I)3-Cl*(U( 1*2668) *(V(I*2686)-V(I*2690)*8.»(V1 1*2689) -V( 1*2607) ) 

2) *V< 1*2688) *(V f 1*2560) -V (1*2816) *8.* (V ( I *2752) -V ( 1*2624) ) ) 

3* W (1*2688 1 MV< 1*1 28) -VI 1*5248) *8 .»( V ( 1*3968) -V (I* 14q8) ) ) 

4*V (1*2686) *U( 1*2686) -V(I*269o)«U(I*269o) *8,4 (V(I*2689)*U( 1+2689) 

5- V(r+2687)»U< 1*258?) ) *V ( I *2560 ) **?.V ( I *2“16 ) *»2+8. * ( V ( J *?752 ) **2 

6- V ( 1*2624) «»f)»V(I+l‘ : 8)4W<I*12S)*V( 1*5248) *W (1*5248) *8,«(V (1*3968) 
7»W(I*3968)-V(I*1408)4H(I* i 4 8 ) ) * V ( 2688 ) **D I V ) 

9*C2*< 16.« (V ( X*2689) +V (I*f6B7)+V( 1*2752) *V (1*2624) *V( 1*3968) *V( I* 
ftl4oB) )-V( 1*2686) -VU*269o)-V( 1*2816) -V( 1*2560 )-V( 1*128) -V( 1*5298) 

B-9o.*V( 1*26881 ) 

1 0 57 0H(J ) =-Cl * (U ( I *268B )* (W { 1*2686) -W( 1*2690) *8, *(K (1*2689) -W( 1*2687) ) 

2) *V ( 1*2688) 4 (W(I*2560)-W( 1*2816) *8.*(W! 1*2752) -W< 1*2624) ) ) 

3*W ( 1*2688) 4 (M ( 1*128) -W ( 1*5248) +&.* (W ( 1 *3968) -W ( 1*1408) ) ) 

4* Wi 1*2686) *U(I*268ft)-W (1*2690) »U( 1*2690) *8,4(«( 1*2689) 4u( 1*2689) 

5-W (1*2687) *U( i* 2687) ) +W( i*2560) *V ( i*2560)-W ( i*2Bl6) *v ( 1*28+ ®> *8. 0 
6(W(I*2752)«V(I*2752)-W(I*‘ i 624)4Vn*2624n*W(I*128)»*2-H(I*5248) 4*2 
7+8, 4 (W( 1*3968) »*2-«( I*l408)**2) +w (2688) 4UIV) 

9*C2M16.4(W( 1*2639) *W (1*2687) *W ( 1*2752) *fl ( I *2624) *W (I *3968) *W ( 1* 

A14P8! ) -W ( 1*2686) -W < I *26?o ) -W ( I*?8' 6) -W< I *2560) -W( 1*128) -W( 1+5248) 

B“90* tt W (1*2688) ) 


1177 


RMS0»RMSD*DIV*42 

1200 


X=U ( 1*2688) *»2 

1202 


RMSU»RMSU*X 

1204 


IFtX.GT.UMAXlUMAXuX 

1213 


Xev ( 1*2888) 4*2 

1215 


RMSV»RMSV*X 

1217 


IF(X.GT.VMAX) VMAX=X 

1223 


X=W(i*2688)442 

1225 


RMSW*RMSW*X 

1227 


IF(X.Gt.WMAX)WMAXpX 

1234 


X=C12*(U(I*2£86) -U(I*269o)*8.«(U(I*2689)-U(I*2687) ) ) 

1243 


SKl=SKl*Xe43 

1246 


SK2=SK2*X**2 

1253 

25 

CONTINUE 

1255 

27 

IEnlE*l 

1257 


IF(M0d(IE. 16).NE.1)50 TO sn 

1263 


DO 30 Ja3 . 1 8 

1264 


DO 30 Ini, 64 

1265 


00 30 L»l*64 

1266 


LL»L*I-1 

1270 


1F(LL.GT.64)LL=LLm64 

1273 

30 

S(I)aS(I)*U(L.Jt3)*U(LL.Jt3) 

1312 

50 

RETURN 

1313 


END 


SUBPROGRAM LENGTH - DVDTMP 
16172 

STATEMENT assignments 

STMT NO* LOCATION STMT NO,* LOCATION STMT NO,* LOCATION STMT NO 

20 ,» 542 27 >* 1256 30 * 12/4 50 

VARIA8LE ASSIGNMENTS 
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START 


The- program START creates' the xnitial velocity field with the 
correct total energy and energy spectrum and a small (but not zero) 
divergence. 
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1 PfjORkAM sTApt (OUT *FsE:t7) 

2 IMPLICIT INTEGER (Z) 

3 lCm/bbi/uh 

A L L M/bN?/(jI 

B DIMENSION VP(2 5b * bA ) .VI (?56,64) 

6 DIMENSION o (f.4,6^,8) , RQl , 2 0 ) ,pQ? (Zq) , 

7 ? VRs(64,64,T , ,VJS ((*4,64,3, . DLG (64 .64 , , 6 , ,Ur ( 64 , j 6 ,64 ) . UI (64 , 

8 3x6,64) ,E(65, ,XK( b5 ) 

9 DIMENSION l7TApE( b ),lTApF(2) 

1 0 EOUIVAlFnCF (VRMI) 'UR(l) )' (VIS'D »UI (1 ))» (f)LG(l) ) 

1) E°t'lVA L ENCF <°U) ,VRn ) ) . < U <163B5> , V I 111 ) 

12 DATA XK/.31 , .44 *.64. .626., 70* .766* .83 • .88 , .941 , .99?, 1.066, 1.175* 

13 21. 2B6, 1.35, 1.4b, I. 55,1. 65, I. 7G, 1 .AS, 2. 0,2. 2, 2, 4, 2. 6 , 2. 8,3.0, 3. 2, 

1 4 33 . 4. 3 . 6 , 3 , 8 , 4. 0*4. 2 . 4. 4, 4. 6 , 4. 8 , 5. 0,5. 5, 6 ,, 6 . 5, 7. ,7. 5, 8 , ,8.5, 9., 

15 410»tll*9l£« *13* » 17« »18*tl9*f?0**30»/ 

lb DATA E/0.f5?»7*l^l*6*101*A*87.5*79*l*67 # 9 # 59*4f5i*^«45*l939.8t 

17 2?9* l*2d#6t 17.7, 13.2 1 lft* • 7 *91 ,6 .32 *5. 06*4 * 09 « 3* 9 2* 07 1 1 . 39, 1 . q 9 , 

1 R 3. 88,. 7c ,.58, . 48 , .j9,.32* ,2 6 * .2?* . l8«.l4, . i i , . 0 RB t , n& 2 . . O 39 , ,0?5, 

19 4, 016, .Oil, . n 08 , .0057, .0042* • 0024,. 0014, .00087, .00053,. 00034. 

20 5.00022, .00014*, 000084.. 000047, .000023, .4 00 01/ 

21 OA T A W1,w2,M/2L*1*2LW2»2LD1'' 

22 DATA R 01 ^ 20 » 0 ,/*RD 2 / 20 # 0 ./ 

23 DATA 17ApE/l ,4 lTARE/,I7TAPE/4,3»0,8LXX0o3593/ 

24 CALL MEMREG (i3l o? 2 * 1 ) 

25 CALL FORqTS(W 1 ,RQ 1 > 

26 CALL FOROTS (W2.RU2) 

27 CA L L CpEATF (Wi.U.RT.O, 0,0, 0*0.0, 1 S 36 ) 

28 CA LL CREATE (W2,U,«T, o,0, 0*0*p, 0*1536) 

29 C=SGRT (3.1415927/ (20.44096. > ) «3 , 1415927e3 4 5flb9 

30 CC= 6 . 28318/20. 

31 RMS=0 • 

32 IX=0 

33 X= 0 . 

34 CALL 0PEN(5lFSET7. 0,2340008. 0,0, 0*3200. tTAPE) 

35 DO 5 0 Z = 1 * 64 

36 X3=7-l 

37 IF (X3.GT.3l.l)x3=X3-64. 

38 DO 2o 3=1 * 64 

39 X2=3-l 

40 IF (X2.GT.3l , 1 ) x2-XH-64. 

4 1 1)0 20 1 = 1 ,64 

42 X1=I-1 

43 IF (Xl.GT.31.DXl=Xl-64. 

44 B=5GR1 (X14a2+X2o42) 

45 IF(R.LT..oo 1 > g O TO 1 ? 

46 p2=XI/8 

47 Pls-X2/B 

48 B=SORT ( (x3«pl ) 442* (X 3 *p 2 ) 4 <f 2 + (X 24 P 1 ) 4 » 2 _ 2 , 4X 1 4 X 2 »Pl *P2 + ( x l»p2 ) 4»2) 

49 G1=-P2«X3/R 

50 02=P1 4X3/P 

51 G3=(P2 tt Xl-Pl«X2)/8 

5? GO TO 13 

53 12 Pl = ). 

54 p2=0, 

55 Gl=0. 

56 Q^sl. 

57 O 3 =0. 

58 13 Y=CC4SQRT(Xl442*x2»42 + x3«42) 

59 DO 1 4 N= 1 *55 

6 q IF(Y.GT.XMn) )G 0 TO j4 
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61 


6? 


63 

14 

64 


66 


66 

15 

67 

16 

68 


69 


7 0 


7] 


72 


73 


74 


76 


76 


77 


r s 


79 


8 o 

20 

81 


8 ? 


83 


84 


86 

3 0 

86 


87 


88 


89 


9 0 

33 

91 


*2 

31 

93 


94 

32 

96 


96 


97 


9R 


99 


100 

5 0 

101 

47 

102 


103 

48 

1 0 A 


105 C 


106 


107 


lo 8 


109 


110 


111 


1 12 


113 


114 

81 

116 


116 

82 

117 


118 

80 

119 


120 



M=N 

60 TO 15 
CONTINUE 
ASO. 

GO TO 16 
A=SORT (E <M> ) 

CONTINUE 

THETA 1=RANF (0. ) *6,2831854 
THeTA2=RAMF (0 .)* 6 .?83l854 

Cl=TOS(THETAl ) 

C2 = GqS (ThEtA?) 

SI =SIN {THETA 1 ) 

S8 =SIN(ThEtAZ> 

G (I *Jtl)=A»(Cl*pI*sl«fll) 

Q (I »J>2)=A»(C 1 «p,J,Sj*02) 

Q ( I , J*3) =**51*03 
0 (I »J»4)s:A*(C2epj +S2*01 ) 

0 (I »Jt5)=Att(C2*P2 + S2*02) 

G (I , J« 6) =A(tS2*Q3 

Cunt inuE 

00 30 K=l,3 

KK=K +3 

DO 30 J=l,64 

CALL EFfg ( 0(1 , J,K) ,Q (, ,J,KK) ,64,) ) 

CONTINUE 
DO 33 K=l,3 
KK=K,3 

00 33 1 = 1,64 

CALL EFT? (0(1*1 ,K) ,G(T,1,KK) ,64,64) 

CONTINUE 
NSEC T= (Z-l ) «24 

1 = IOONE (V)I ) 

IF (I. NE. DRO TO 3| 

I=IOONE (k2j 

IFM.NE.DGO To 32 

SMALL 0UT(Q(1 »1 ,1) *DLG(1,1 *4) ,12288) 

SMALL OUT (0(1 *1 « 41 *°LP(1 ,1 » 7 ) .J2288) 

I = IRANW(RQ 1 .DLG(L,JL,4) , 122 H 8,NSECT,W) ) 

1= Ip A nIk (RQ2,DL6(I»I,7) ,122 ea ,NSECT,W2) 

CONTINUE 
I=I00HE(wl ) 

IF (I ,NE.l )G0 To at 
I=ID 0NE(w2) 

IF ( I ,NE, ) ) 60 To 48 

READ 6A j s BY I& J s BY 64 z s INTO L a oGE CORE 
NINC= (I 

DO 200 N=l,4 
NSECT)=NINC 
DO lbO M-1,3 
NSEC T-NSECT ) 

DO 80 Z-l ,64 

I=IrANR(RQj,uR(i, 1 .Z> ,i024,N$ECT,W 1 ) 

I=IR A nR (RQs.UI ( 1,1, Z) ,1024, NSFCT.W2) 

1= IDONE ( WI ) 

IF (I ,NE.l ! GD To 8l 
I=ID0NE(W2) 

IF ( I ,NE, 1 ) GO To 82 

NSECT=NSECT,24 

DO 95 J=l,l6,4 
DO 9 0 Z=j ,64 
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121 

12? 

90 

123 

124 

91 

125 

]2& 

12 7 

^2 

128 

95 

129 

13n 

131 

132 

101 

133 

134 

IOO 

135 

150 

136 

2u0 

137 
1 3fi 
139 
1*0 

321 

1 A 1 

142 

143 

144 
!45 
146 
1*7 

3()0 

1*8 

3 1 0 

1-+9 

320 

1 bO 

151 

152 

350 

153 

154 
lb5 



Sf-’ALL IN (VP(1,2) ,UR(1 ,J,Z) ,2*6) 
SMALL IN (V T <l,Z) .UI (] ,J,Z) ,?56) 

DO 91 1=1, ?C6 

CALL FF\Z ( VP ( I * 1 ) * V I ( I , 1 ) ,64,256) 

00 9 2 Z=j ,64 

SMALL OUT ( Vr ( 1 » Z ) * UR ( 1 * J , Z) ,p56 ) 
SMALL 0UT1VT (lfZ) »UI (, , J,Z) » ? 56) 
CONTINUE 
NSecT=NSecT) 

DO 109 Z=I,6* 

1 = IRA:MW (RQJ ,uR(l,i,Z) ,192*»NSECT.W1 ) 
I*I00nE(W1) 

IF(I.Nt.l)GO To 101 
NSECT=NSECT+Z-» 

NSECTl=NSECTl+8 

NINC=«iINC+ 2 
DO 320 Ksl*1S 
NSErT=(K-n«i2e 

I=IRA9K(RQi.uRf 1 ) , 65536, nSECt,W 1) 

I = ID0ME ( Mil ) 

IF<I.'JE.1)S 0 To 321 

DO 3iO H=l,? 

Il=32 /68 <MM_ilu 

SMALL IN(0(D ,UR!II> ,32768) 

DO 3 0 , I=i, 32 7 68 
Q 1 1 ) =C*t) ( I ) 

RNS=RmS+Q(I)-»»2 
WfirTE(7) , (0(1) ,1=1,32768, 

CONTINUE 

RmS=S0RT,R m s/3.)/51?. 
pRINT 3bO , RMS 
FORMAT („ RMS=#E14,7) 

CALL AFSREl (5 LFSETT,o.i7tAPE) 

STOP 

ENO 
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